class: center, middle, inverse, title-slide # Displaying analytic results ### Abhijit Dasgupta ### BIOF 339 --- class: middle,center,inverse # Comparing groups --- layout: true ## The `ggpubr` package --- The **ggpubr** package, which extends **ggplot2** functionality, has several functions that allow the computation and visualization of different statistical analyses Under the hood, it's just fancy application of R for statistical tests and then translating the results to ggplot geometries. --- .pull-left[ ```r plt <- ggboxplot(penguins, x = 'species', y = 'body_mass_g', color = 'species', add='jitter', add.params = list(color='grey', size=0.5)) summary.stats <- penguins %>% select(body_mass_g, species) %>% group_by(species) %>% get_summary_stats(type='common') plt2 <- ggsummarytable( summary.stats, x = 'species', y = c('n','median','iqr') ) + theme_minimal()+ theme(panel.grid=element_blank(), axis.text.x = element_blank())+ labs(x='', y='') ggarrange(plt, plt2, ncol=1, heights=c(3,1)) ``` ] .pull-right[ <!-- --> ] --- .pull-left[ ```r (plt <- ggplot(penguins, aes(x=species, y=body_mass_g, color=species))+ geom_boxplot()+ geom_jitter(color='grey', size=0.5)+ theme_classic()+ theme(legend.position = 'top')) ``` ] .pull-right[ <!-- --> ] --- .pull-left[ ```r (plt <- ggplot(penguins, aes(x=species, y=body_mass_g, color=species))+ geom_boxplot()+ geom_jitter(color='grey', size=0.5)+ theme_classic()+ theme(legend.position = 'top')) ``` ```r summary.stats <- penguins %>% select(body_mass_g, species) %>% group_by(species) %>% get_summary_stats(type='common') (plt2 <- ggsummarytable( summary.stats, x = 'species', y = c('n','median','iqr'), color='species' ) + theme_minimal()+ theme(panel.grid=element_blank(), axis.text.x = element_blank())+ labs(x='', y='')) ``` ] .pull-right[ <!-- --> ] --- .pull-left[ ```r (plt <- ggplot(penguins, aes(x=species, y=body_mass_g, color=species))+ geom_boxplot()+ geom_jitter(color='grey', size=0.5)+ theme_classic()+ theme(legend.position = 'top')) ``` ```r summary.stats <- penguins %>% select(body_mass_g, species) %>% group_by(species) %>% get_summary_stats(type='common') (plt2 <- ggsummarytable( summary.stats, x = 'species', y = c('n','median','iqr'), color='species' ) + theme_minimal()+ theme(panel.grid=element_blank(), axis.text.x = element_blank())+ labs(x='', y='')) ``` ] .pull-right[ ```r ggarrange(plt, plt2, ncol=1, heights = c(4,1)) ``` <!-- --> ] --- .pull-left[ ```r (summ_plt <- ggsummarystats( penguins, x = 'species', y = 'body_mass_g', ggfunc = ggboxplot, add='jitter', color='species', add.params=list(color='grey', size=0.5) )) ``` ] .pull-right[ <!-- --> ] --- .pull-left[ ```r (summ_plt <- ggsummarystats( penguins, x = 'species', y = 'body_mass_g', ggfunc = ggboxplot, add='jitter', color='species', add.params=list(color='grey', size=0.5) )) ``` ```r my_comparisons <- list(c('Adelie','Chinstrap'), c('Chinstrap','Gentoo'), c('Adelie','Gentoo')) summ_plt$main.plot <- summ_plt$main.plot + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y = 8000) summ_plt ``` ----- You can also add facets using the argument `facet.by` in `ggsummarystats` ] .pull-right[ <!-- --> ] --- count: false Scatter plots .panel1-corr-auto[ ```r *ggscatter(penguins, * x = 'bill_length_mm', * y = 'body_mass_g', * color='grey', * add='reg.line', * add.params=list(color='blue')) ``` ] .panel2-corr-auto[ <!-- --> ] --- count: false Scatter plots .panel1-corr-auto[ ```r ggscatter(penguins, x = 'bill_length_mm', y = 'body_mass_g', color='grey', add='reg.line', add.params=list(color='blue')) + * labs(x = 'Bill length (mm)', * y = 'Body mass (g)') ``` ] .panel2-corr-auto[ <!-- --> ] --- count: false Scatter plots .panel1-corr-auto[ ```r ggscatter(penguins, x = 'bill_length_mm', y = 'body_mass_g', color='grey', add='reg.line', add.params=list(color='blue')) + labs(x = 'Bill length (mm)', y = 'Body mass (g)') + * stat_cor( label.x=30, label.y=6000, * label.sep='\n') ``` ] .panel2-corr-auto[ <!-- --> ] --- count: false Scatter plots .panel1-corr-auto[ ```r ggscatter(penguins, x = 'bill_length_mm', y = 'body_mass_g', color='grey', add='reg.line', add.params=list(color='blue')) + labs(x = 'Bill length (mm)', y = 'Body mass (g)') + stat_cor( label.x=30, label.y=6000, label.sep='\n') + * facet_wrap(~species) ``` ] .panel2-corr-auto[ <!-- --> ] <style> .panel1-corr-auto { color: black; width: 38.6060606060606%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel2-corr-auto { color: black; width: 59.3939393939394%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel3-corr-auto { color: black; width: NA%; hight: 33%; float: left; padding-left: 1%; font-size: 80% } </style> --- ### Add tables to a graphic .pull-left[ ```r library(ggsci) # Themes for science journals dens_plt <- ggplot(penguins, aes(x = body_mass_g))+ geom_density(aes(fill=species))+ scale_fill_jco(alpha=0.3)+ # Journal of Clinical Oncology colors theme_classic() + theme(axis.text.y = element_blank(), legend.position='top') + labs(y="") stable <- desc_statby(penguins, measure.var='body_mass_g', grps='species') stable <- stable[,c('species','length','mean','sd')] stable_plt <- ggtexttable(stable, rows=NULL) ggarrange(dens_plt, stable_plt, ncol=1, heights=c(2,1)) ``` ] .pull-right[ <!-- --> ] --- ### Add tables to a graphic .pull-left[ ```r library(ggsci) dens_plt <- ggplot(penguins, aes(x = body_mass_g))+ geom_density(aes(fill=species))+ scale_fill_jco(alpha=0.3)+ theme_classic() + theme(axis.text.y = element_blank(), legend.position='top') + * stat_central_tendency( * aes(color=species), * type='mean', * geom='line', linetype=2 * ) + * scale_color_jco() + labs(y = "") stable <- desc_statby(penguins, measure.var='body_mass_g', grps='species') stable <- stable[,c('species','length','mean','sd')] stable_plt <- ggtexttable(stable, rows=NULL) ggarrange(dens_plt, stable_plt, ncol=1, heights=c(2,1)) ``` ] .pull-right[ <!-- --> ] --- .pull-left[ ```r ggscatter(penguins, x='bill_length_mm', y='body_mass_g', color='species', shape='species', size=0.5, ellipse=TRUE)+ stat_mean(aes(color=species, shape=species), size=4) ``` ] .pull-right[ <!-- --> ] --- layout: true ## Modeling results --- Let's start with a fairly naive multiple regression model ```r m1 <- lm(body_mass_g ~., data=penguins) summary(m1)$coef ``` ``` Estimate Std. Error t value Pr(>|t|) (Intercept) 84087.94460 41912.018942 2.0062967 4.565793e-02 speciesChinstrap -282.53947 88.790344 -3.1820968 1.604233e-03 speciesGentoo 890.95818 144.563079 6.1631102 2.124019e-09 islandDream -21.18017 58.390407 -0.3627337 7.170410e-01 islandTorgersen -58.77708 60.852126 -0.9659003 3.348167e-01 bill_length_mm 18.96388 7.111843 2.6665210 8.049880e-03 bill_depth_mm 60.79775 20.002232 3.0395481 2.562907e-03 flipper_length_mm 18.50387 3.128430 5.9147481 8.464766e-09 sexmale 378.97704 48.074295 7.8831534 4.948533e-14 year -42.78456 20.949432 -2.0422775 4.193578e-02 ``` --- count: false .panel1-mod1-auto[ ```r *library(gtsummary) ``` ] .panel2-mod1-auto[ ] --- count: false .panel1-mod1-auto[ ```r library(gtsummary) *theme_gtsummary_compact() ``` ] .panel2-mod1-auto[ ] --- count: false .panel1-mod1-auto[ ```r library(gtsummary) theme_gtsummary_compact() *tbl_regression(m1, * label = list( * bill_length_mm ~ 'Bill length (mm)', * bill_depth_mm ~ 'Bill depth (mm)', * flipper_length_mm ~ 'Flipper length (mm)', * species ~ 'Species', * island ~ 'Island', * sex ~ 'Sex', * year ~ 'Year') *) ``` ] .panel2-mod1-auto[ preservedb302acb7dea8a3d ] --- count: false .panel1-mod1-auto[ ```r library(gtsummary) theme_gtsummary_compact() tbl_regression(m1, label = list( bill_length_mm ~ 'Bill length (mm)', bill_depth_mm ~ 'Bill depth (mm)', flipper_length_mm ~ 'Flipper length (mm)', species ~ 'Species', island ~ 'Island', sex ~ 'Sex', year ~ 'Year') ) %>% * add_global_p() ``` ] .panel2-mod1-auto[ preserve70e3efed31ef5a1c ] --- count: false .panel1-mod1-auto[ ```r library(gtsummary) theme_gtsummary_compact() tbl_regression(m1, label = list( bill_length_mm ~ 'Bill length (mm)', bill_depth_mm ~ 'Bill depth (mm)', flipper_length_mm ~ 'Flipper length (mm)', species ~ 'Species', island ~ 'Island', sex ~ 'Sex', year ~ 'Year') ) %>% add_global_p() %>% * bold_p(t = 0.05) ``` ] .panel2-mod1-auto[ preservee538e67da83dc682 ] --- count: false .panel1-mod1-auto[ ```r library(gtsummary) theme_gtsummary_compact() tbl_regression(m1, label = list( bill_length_mm ~ 'Bill length (mm)', bill_depth_mm ~ 'Bill depth (mm)', flipper_length_mm ~ 'Flipper length (mm)', species ~ 'Species', island ~ 'Island', sex ~ 'Sex', year ~ 'Year') ) %>% add_global_p() %>% bold_p(t = 0.05) %>% * bold_labels() ``` ] .panel2-mod1-auto[ preserve4aae10d824418213 ] --- count: false .panel1-mod1-auto[ ```r library(gtsummary) theme_gtsummary_compact() tbl_regression(m1, label = list( bill_length_mm ~ 'Bill length (mm)', bill_depth_mm ~ 'Bill depth (mm)', flipper_length_mm ~ 'Flipper length (mm)', species ~ 'Species', island ~ 'Island', sex ~ 'Sex', year ~ 'Year') ) %>% add_global_p() %>% bold_p(t = 0.05) %>% bold_labels() %>% * italicize_levels() ``` ] .panel2-mod1-auto[ preserve3f46698eb29c82dc ] <style> .panel1-mod1-auto { color: black; width: 38.6060606060606%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel2-mod1-auto { color: black; width: 59.3939393939394%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel3-mod1-auto { color: black; width: NA%; hight: 33%; float: left; padding-left: 1%; font-size: 80% } </style> --- ### Putting together multiple models .pull-left[ ```r (gt_r1 <- glm(response ~ trt + grade, trial, family = binomial) %>% tbl_regression(exponentiate = TRUE)) ``` ] .pull-right[ preservece78e1b69e63de5b ] --- ### Putting together multiple models .pull-left[ ```r (gt_r2 <- coxph(Surv(ttdeath, death) ~ trt + grade, trial) %>% tbl_regression(exponentiate = TRUE)) ``` ] .pull-right[ preservea8e8d46170a799d7 ] --- ### Putting together multiple models .pull-left[ ```r (gt_t1 <- trial[c("trt", "grade")] %>% tbl_summary(missing = "no") %>% add_n() %>% modify_header(stat_0 ~ "**n (%)**") %>% modify_footnote(stat_0 ~ NA_character_)) ``` ] .pull-right[ preserve3e4bd5fb1314577c ] --- ### Putting together multiple models ```r tbl_merge( list(gt_t1, gt_r1, gt_r2), tab_spanner = c(NA_character_, "**Tumor Response**", "**Time to Death**") ) ``` preservecbab918c1ae4dd6d --- ### The finalfit package ```r explanatory <- setdiff(names(penguins), 'body_mass_g') dependent <- 'body_mass_g' t2 <- penguins %>% finalfit::finalfit(dependent, explanatory, metrics=TRUE) knitr::kable(t2[[1]], row.names = F) ``` <table> <thead> <tr> <th style="text-align:left;"> Dependent: body_mass_g </th> <th style="text-align:left;"> </th> <th style="text-align:left;"> unit </th> <th style="text-align:left;"> value </th> <th style="text-align:left;"> Coefficient (univariable) </th> <th style="text-align:left;"> Coefficient (multivariable) </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> species </td> <td style="text-align:left;"> Adelie </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 3700.7 (458.6) </td> <td style="text-align:left;"> - </td> <td style="text-align:left;"> - </td> </tr> <tr> <td style="text-align:left;"> </td> <td style="text-align:left;"> Chinstrap </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 3733.1 (384.3) </td> <td style="text-align:left;"> 32.43 (-100.37 to 165.22, p=0.631) </td> <td style="text-align:left;"> -282.54 (-457.22 to -107.86, p=0.002) </td> </tr> <tr> <td style="text-align:left;"> </td> <td style="text-align:left;"> Gentoo </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 5076.0 (504.1) </td> <td style="text-align:left;"> 1375.35 (1264.91 to 1485.80, p<0.001) </td> <td style="text-align:left;"> 890.96 (606.55 to 1175.36, p<0.001) </td> </tr> <tr> <td style="text-align:left;"> island </td> <td style="text-align:left;"> Biscoe </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4716.0 (782.9) </td> <td style="text-align:left;"> - </td> <td style="text-align:left;"> - </td> </tr> <tr> <td style="text-align:left;"> </td> <td style="text-align:left;"> Dream </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 3712.9 (416.6) </td> <td style="text-align:left;"> -1003.11 (-1149.16 to -857.07, p<0.001) </td> <td style="text-align:left;"> -21.18 (-136.05 to 93.69, p=0.717) </td> </tr> <tr> <td style="text-align:left;"> </td> <td style="text-align:left;"> Torgersen </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 3706.4 (445.1) </td> <td style="text-align:left;"> -1009.65 (-1206.75 to -812.54, p<0.001) </td> <td style="text-align:left;"> -58.78 (-178.49 to 60.94, p=0.335) </td> </tr> <tr> <td style="text-align:left;"> bill_length_mm </td> <td style="text-align:left;"> [32.1,59.6] </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4201.8 (802.0) </td> <td style="text-align:left;"> 87.42 (74.82 to 100.01, p<0.001) </td> <td style="text-align:left;"> 18.96 (4.97 to 32.96, p=0.008) </td> </tr> <tr> <td style="text-align:left;"> bill_depth_mm </td> <td style="text-align:left;"> [13.1,21.5] </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4201.8 (802.0) </td> <td style="text-align:left;"> -191.64 (-229.84 to -153.45, p<0.001) </td> <td style="text-align:left;"> 60.80 (21.45 to 100.15, p=0.003) </td> </tr> <tr> <td style="text-align:left;"> flipper_length_mm </td> <td style="text-align:left;"> [172.0,231.0] </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4201.8 (802.0) </td> <td style="text-align:left;"> 49.69 (46.70 to 52.67, p<0.001) </td> <td style="text-align:left;"> 18.50 (12.35 to 24.66, p<0.001) </td> </tr> <tr> <td style="text-align:left;"> sex </td> <td style="text-align:left;"> female </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 3862.3 (666.2) </td> <td style="text-align:left;"> - </td> <td style="text-align:left;"> - </td> </tr> <tr> <td style="text-align:left;"> </td> <td style="text-align:left;"> male </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4545.7 (787.6) </td> <td style="text-align:left;"> 683.41 (526.02 to 840.80, p<0.001) </td> <td style="text-align:left;"> 378.98 (284.40 to 473.56, p<0.001) </td> </tr> <tr> <td style="text-align:left;"> year </td> <td style="text-align:left;"> 2007 </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4124.5 (792.5) </td> <td style="text-align:left;"> - </td> <td style="text-align:left;"> - </td> </tr> <tr> <td style="text-align:left;"> </td> <td style="text-align:left;"> 2008 </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4266.7 (789.5) </td> <td style="text-align:left;"> - </td> <td style="text-align:left;"> - </td> </tr> <tr> <td style="text-align:left;"> </td> <td style="text-align:left;"> 2009 </td> <td style="text-align:left;"> Mean (sd) </td> <td style="text-align:left;"> 4210.3 (822.9) </td> <td style="text-align:left;"> - </td> <td style="text-align:left;"> - </td> </tr> <tr> <td style="text-align:left;"> NA </td> <td style="text-align:left;"> NA </td> <td style="text-align:left;"> NA </td> <td style="text-align:left;"> NA </td> <td style="text-align:left;"> 41.42 (-63.17 to 146.02, p=0.437) </td> <td style="text-align:left;"> -42.78 (-84.00 to -1.57, p=0.042) </td> </tr> </tbody> </table> --- ### Other packages + **sjPlot** [(link)](https://strengejacke.github.io/sjPlot/articles/tab_model_estimates.html) + **stargazer** #### Utility packages The function `broom::tidy` takes results from [many](https://broom.tidymodels.org/articles/available-methods.html) models and transforms them into _tidy datasets_, which can be formatted and output using a variety of other packages in the R ecosystem. Some popular and effective packages are + **gt**: The grammar of tables [(link)](https://gt.rstudio.com) + **flextable**: Creating tables that can work with RMarkdown/HTML as well as the [officeverse](https://ardata-fr.github.io/officeverse/), which creates Microsoft Word and Powerpoint files from R and RMarkdown [(link)](https://ardata-fr.github.io/flextable-book/index.html) + **pixiedust** is designed to format the results from the `broom::tidy` [(link)](https://github.com/nutterb/pixiedust) * **xtable** and **stargazer** are targeted towards \LaTeX and PDF outputs .footnote[This slide differs from the video. It was updated in 2021] --- layout: true ## Model results ### Plotting results --- .pull-left[ ```r out <- broom::tidy(m1) %>% slice(-1) %>% mutate(lcb = estimate - 2*std.error, ucb = estimate + 2*std.error) %>% select(term, estimate, lcb, ucb) knitr::kable(out, digits = 2) %>% kable_styling() ``` ] .pull-right[ <table class="table" style="margin-left: auto; margin-right: auto;"> <thead> <tr> <th style="text-align:left;"> term </th> <th style="text-align:right;"> estimate </th> <th style="text-align:right;"> lcb </th> <th style="text-align:right;"> ucb </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> speciesChinstrap </td> <td style="text-align:right;"> -282.54 </td> <td style="text-align:right;"> -460.12 </td> <td style="text-align:right;"> -104.96 </td> </tr> <tr> <td style="text-align:left;"> speciesGentoo </td> <td style="text-align:right;"> 890.96 </td> <td style="text-align:right;"> 601.83 </td> <td style="text-align:right;"> 1180.08 </td> </tr> <tr> <td style="text-align:left;"> islandDream </td> <td style="text-align:right;"> -21.18 </td> <td style="text-align:right;"> -137.96 </td> <td style="text-align:right;"> 95.60 </td> </tr> <tr> <td style="text-align:left;"> islandTorgersen </td> <td style="text-align:right;"> -58.78 </td> <td style="text-align:right;"> -180.48 </td> <td style="text-align:right;"> 62.93 </td> </tr> <tr> <td style="text-align:left;"> bill_length_mm </td> <td style="text-align:right;"> 18.96 </td> <td style="text-align:right;"> 4.74 </td> <td style="text-align:right;"> 33.19 </td> </tr> <tr> <td style="text-align:left;"> bill_depth_mm </td> <td style="text-align:right;"> 60.80 </td> <td style="text-align:right;"> 20.79 </td> <td style="text-align:right;"> 100.80 </td> </tr> <tr> <td style="text-align:left;"> flipper_length_mm </td> <td style="text-align:right;"> 18.50 </td> <td style="text-align:right;"> 12.25 </td> <td style="text-align:right;"> 24.76 </td> </tr> <tr> <td style="text-align:left;"> sexmale </td> <td style="text-align:right;"> 378.98 </td> <td style="text-align:right;"> 282.83 </td> <td style="text-align:right;"> 475.13 </td> </tr> <tr> <td style="text-align:left;"> year </td> <td style="text-align:right;"> -42.78 </td> <td style="text-align:right;"> -84.68 </td> <td style="text-align:right;"> -0.89 </td> </tr> </tbody> </table> ] --- .pull-left[ ```r (plt1 <- out %>% ggplot(aes(x = term, y = estimate, ymin = lcb, ymax = ucb))+ geom_pointrange()+ geom_hline(yintercept=0, linetype=2)+ xlab('') + geom_vline(xintercept=0, linetype=2)+ coord_flip() + theme_classic()) ``` ] .pull-right[ <!-- --> ] --- .pull-left[ ```r (plt1 <- out %>% ggplot(aes(x = term, y = estimate, ymin = lcb, ymax = ucb))+ geom_pointrange()+ geom_hline(yintercept=0, linetype=2)+ xlab('') + geom_vline(xintercept=0, linetype=2)+ coord_flip() + theme_classic()) ``` ```r out <- out %>% mutate(across(-term, round, 2)) %>% mutate(ci = glue::glue('{estimate} ({lcb},{ucb})')) plt2 <- ggplot(out, aes(x = term, y = 0))+ geom_text(aes(label=ci), size=3, hjust=0)+ coord_flip()+theme_void() + scale_y_continuous(limits = c(-0.1, 0.2)) ggpubr::ggarrange(plt1, plt2, nrow=1, widths=c(2,1), align='h') ``` ] .pull-right[ <!-- --> ] --- We can use the **coefplot** package to make a similar plot .pull-left[ ```r pacman::p_load("coefplot") coefplot(m1, intercept=FALSE) ``` .footnote[This slide differs from the video. It was updated in 2021] ] .pull-right[ <!-- --> ] --- .pull-left[ We can also use the **see** package ```r pacman::p_load(see, parameters) plot(model_parameters(m1), show_labels = TRUE) ``` .footnote[This slide differs from the video. It was updated in 2021] ] .pull-right[ <!-- --> ] --- layout: true ## Model results ### Survival analysis --- .pull-left[ ```r library(survival) library(survminer) pbc <- pbc %>% mutate(trt = as.factor(trt), stage = as.factor(stage)) m <- survfit(Surv(time, status==2)~trt, data=pbc) ggsurvplot(m, pval=TRUE, risk.table=T) ``` ] .pull-right[ <!-- --> ] --- .pull-left[ ```r library(survival) library(survminer) pbc <- pbc %>% mutate(trt = as.factor(trt), stage = as.factor(stage)) m <- survfit(Surv(time, status==2)~trt, data=pbc) ggsurvplot(m, pval=TRUE, risk.table=T) ``` ```r m2 <- coxph(Surv(time, status==2) ~ trt, data = pbc) gtsummary::tbl_regression(m2, exponentiate=T) ``` ] .pull-right[ preservefa538ab9eac99275 ] --- .pull-left[ ```r m3 <- coxph(Surv(time, status==2) ~ trt + sex + age + stage, data=pbc) out <- broom::tidy(m3) %>% mutate(lcb = estimate - 2*std.error, ucb = estimate + 2*std.error) %>% mutate(across(c(estimate,lcb,ucb), exp)) *gtsummary::tbl_regression(m3, exponentiate = T) ``` preserve5d634018ef9b25bf ] .pull-right[ ```r out %>% ggplot(aes(x = term, y = estimate, ymin = lcb, ymax=ucb))+ geom_pointrange()+ geom_hline(yintercept=1, linetype=2)+ coord_flip()+ theme_classic() ``` <!-- --> ] --- .pull-left[ ```r ggsurvplot( m, pval = TRUE, risk.table = FALSE, palette="jama", xlab = "Days", legend.title="", * facet.by = "stage", legend.labs=c("Control","Treatment") ) ``` .footnote[This slide differs from the video and was updated in 2021] ] .pull-right[ <!-- --> ]