Different ways to estimate treatment effects

in longitudinal randomised controlled trials

Agustin Calatroni
AXC

Tuesday, the 19th of March, 2024

Abstract

Regarding the analysis of longitudinal RCT data there is a debate going on whether an adjustment for the baseline value of the outcome variable should be made. When an adjustment is made, there is a lot of misunderstanding regarding the way this should be done. Therefore, the aims to:

  • to explain different methods used to estimate treatment effects in RCTs

  • to illustrate the different methods with a real data example

  • to give an advise on how to analyse RCT data

Example dataset

Effectiveness of a long-term homocysteine-lowering treatment with folic acid plus pyridoxine in reducing systolic blood.

Tables

In this 24 months , randomized, placebo-controlled trial, a baseline measurement and two follow-up measurements (after 12 and after 24 months) were performed. At each time-point systolic blood pressure was measured four times and the average value was used in the analysis.

Code 
add_by_n <- function(data, variable, by, tbl, ...) {
   data %>%
      select(all_of(c(variable, by))) %>%
      dplyr::group_by(.data[[by]]) %>%
      dplyr::summarise_all(~sum(!is.na(.))) %>%
      rlang::set_names(c("by", "variable")) %>%
      dplyr::left_join(
         tbl$df_by %>% select(by, by_col),
         by = "by"
      ) %>%
      mutate(
         by_col = paste0("add_n_", by_col),
         variable = style_number(variable)
      ) %>%
      select(-by) %>%
      tidyr::pivot_wider(names_from = by_col, 
                         values_from = variable)}

tbl_summary(data = dat_twisk_w %>% mutate(TRT01P = fct_rev(TRT01P)),
            by = TRT01P , 
            include = starts_with('SYS_'),
            digit = list(everything() ~ c(1,1)),
            statistic = list(all_continuous() ~ "{mean} ({sd})"),
            missing = 'no',
            label = list( SYS_0 ~ 'Week 0 (Randomization)',
                          SYS_1 ~ 'Week 12',
                          SYS_2 ~ 'Week 24')) %>% 
   add_stat( fns = everything() ~ add_by_n ) %>%
   add_difference(test = list(all_continuous() ~ "t.test"),
                  pvalue_fun = ~style_pvalue(.x, digits = 2),
                  estimate_fun = list(all_continuous() ~ function(x) style_number(x, digits=1)) ) %>% 
   bold_labels() %>% 
   modify_column_merge(pattern = "{estimate} ({ci})") %>% 
   modify_header(
      label = '',
      estimate = '**DIFF (95% CI)**',
      all_stat_cols(FALSE) ~ "**{level}**<br>n = {n}",
      starts_with("add_n_stat") ~ "**n**") %>% 
   modify_table_body(
      ~ .x %>%
         dplyr::relocate(add_n_stat_1, .before = stat_1) %>%
         dplyr::relocate(add_n_stat_2, .before = stat_2)
   ) %>% 
   modify_spanning_header(
      c(all_stat_cols(F), starts_with("add_n_stat")) ~ "**Treatment**") %>% 
   as_gt() %>% 
   gtsave('tbl/tbl_1.png', zoom = 4)

knitr::include_graphics('tbl/tbl_1.png')

Code 
pacman::p_load(gtreg)

tbl_reg_summary(data = dat_twisk_w %>% mutate(TRT01P = fct_rev(TRT01P)),
                by = TRT01P , 
                include = starts_with('SYS_'),
                label = list( SYS_0 ~ 'Week 0 (Randomization)',
                              SYS_1 ~ 'Week 12',
                              SYS_2 ~ 'Week 24'),
                digits = everything() ~ c(0, 1, 1, 0, 0, 0, 0, 0) ) %>% 
   add_overall() %>% 
   bold_labels() %>% 
   modify_header(
      label = '',
      all_stat_cols(TRUE) ~ "**{level}**" ) %>% 
   modify_spanning_header(
      c(all_stat_cols(F), starts_with("add_n_stat")) ~ "**Treatment**")%>% 
   as_gt() %>% 
   gtsave('tbl/tbl_2.png', zoom = 3)

knitr::include_graphics('tbl/tbl_2.png')

Figures

Code 
ggplot(data = dat_twisk_l %>% drop_na(AVAL),
       aes(x = AVISITN, y = AVAL, color = TRT01P)) +
   geom_quasirandom(dodge.width = 0.5, alpha = 0.5, size = 2) +
   scale_colour_brewer(palette = "Set1",
                       name = 'Treatment') +
   scale_y_continuous(name = 'Systolic Blood Pressure (mmHg)',
                      limits = c(80, 200)) +
   scale_x_continuous(name = '',
                      breaks = c(0, 1, 2),
                      labels = c('Week 0\nRandomization', 'Week 12', 'Week 24')) +
   theme_bw(base_size = 18) + 
   theme(panel.grid.minor = element_blank(),
         legend.position      = c(0.99, 0.99),
         legend.justification = c('right','top'))

Code 
m <- mmrm(AVAL ~ TRT01P * AVISIT + us(AVISIT | USUBJID), data = dat_twisk_l)

e <- emmeans(m, revpairwise ~ TRT01P | AVISIT, adjust = 'none') %>% 
   as.data.frame() %>% 
   filter(contrast != '.') %>% 
   data.frame()

f_a <- ggplot(data = dat_twisk_l %>% drop_na(AVAL),
              aes(x = AVISITN, y = AVAL, color = TRT01P)) +
   geom_quasirandom(dodge.width=1, alpha = 0.2, size = 2.5, show.legend = FALSE) +
   geom_boxplot(aes(group = interaction(AVISITN, TRT01P)),
                fill = 'transparent',
                outlier.shape = NA,
                show.legend = FALSE) +
   scale_colour_brewer(palette = "Set1",
                       name = 'Treatment') +
   stat_summary(fun = mean,
                geom = "point",
                size = 3,
                position = position_dodge(width = 1)) +
   scale_y_continuous(name = 'Systolic Blood Pressure (mmHg)',
                      limits = c(80, 200)) +
   scale_x_continuous(name = NULL,
                      breaks = c(0, 1, 2),
                      labels = c('Week 0 (R)', 'Week 12', 'Week 24')) +
   theme_bw(base_size = 14) + 
   theme(panel.grid.minor = element_blank(),
         legend.position      = c(0.99, 0.99),
         legend.justification = c('right','top'))

f_b <- ggplot(data = e,
              aes(x = AVISIT, y = emmean, group = contrast)) + 
   geom_hline(yintercept  = 0, color = 'gray50') +
   geom_point(size = 4) +
   geom_line(linewidth = 1) +
   geom_linerange(aes(ymin = lower.CL,
                      ymax = upper.CL),
                  size = 1.5) +
   geom_text(aes( y = -14.5,
                  label = gtsummary::style_pvalue(p.value, prepend_p = TRUE) ) ) +
   scale_color_discrete(name = 'Treatment') +
   scale_y_continuous(name = 'DIFF (95% CI)',
                      limits = c(-15, 2)) +
   scale_x_discrete(name = '') +
   theme_bw(base_size = 14) +
   theme(panel.grid.minor = element_blank(),
         axis.text.x=element_blank())

f_a / f_b + plot_layout( heights = c(2, 1))

Code 
ggplot(data = dat_twisk_l %>% drop_na(AVAL),
       aes(x = AVISIT, y = AVAL, color = TRT01P , group = USUBJID )) +
   geom_quasirandom(dodge.width = 0.5, alpha = 0.5, size = 2.5) +
   geom_line(alpha = 0.5,
             position=position_quasirandom(dodge.width = 0.5)) +
   facet_wrap(. ~ TRT01P  ) +
   scale_colour_brewer(palette = "Set1",
                       name = 'Treatment') +
   scale_y_continuous(name = 'Systolic Blood Pressure (mmHg)',
                      limits = c(80, 200)) +
   labs(x = NULL) +
   theme_bw(base_size = 18) + 
   theme(panel.grid.minor = element_blank(),
         legend.position      = c(0.99, 0.99),
         legend.justification = c('right','top'))

trelliscope: Create and Explore Data Frames of Visualizations

Trelliscope is an open source project that provides tools for data scientists to build and share interactive collections of visualizations for detailed data exploration and more.

R/trelliscope Code 
pacman::p_load(trelliscope)

summary_df <- dat_twisk_l %>% 
   select(USUBJID, TRT01P, AVISITN, AVISIT, AVAL, CHG) %>%
   drop_na(AVAL) %>% 
   summarise(
      mean_AVAL = mean(AVAL),
      sd_AVAL = sd(AVAL),
      max_AVAL  = max(AVAL),
      mean_CHG  = mean(CHG),
      sd_CHG  = sd(CHG),
      max_CHG  = max(CHG),
      .by = c(USUBJID, TRT01P))

panels_df <- (
   ggplot(dat_twisk_l %>% drop_na(AVAL) %>% mutate(TRT01P = factor(TRT01P, levels = c('TRT', 'CNT'))),
          aes(x = AVISIT, y = AVAL, color = TRT01P, group = USUBJID) ) +
      geom_point(size = 4) +
      geom_line(linewidth = 1) +
      scale_colour_manual(values = c( 'TRT'  = '#377eb8',
                                      'CNT'  = '#e41a1c') ) +
      scale_y_continuous(name = 'Systolic Blood Pressure (mmHg)',
                         limits = c(80, 200)) +
      scale_x_discrete(name = NULL) +
      theme_bw(base_size = 14) +
      theme(legend.position = "none",
            panel.grid.minor = element_blank()) +
      facet_panels( vars(USUBJID, TRT01P), scales = 'same') ) %>% 
   as_panels_df()

trellis_df <- left_join(summary_df, 
                        panels_df, 
                        by = c("USUBJID", "TRT01P")) %>% 
   as_trelliscope_df(
      name = "Explore Individual Subject Data",
      path = './_trelliscope',
      description = "Effectiveness of a long-term homocysteine-lowering treatment 
      with folic acid plus pyridoxine in reducing systolic blood")%>%
   set_default_labels(c("USUBJID", "TRT01P")) %>% 
   set_default_layout(ncol = 3, page = 1) %>% 
   set_var_labels(
      mean_AVAL = 'Mean SBP') %>% 
   set_default_sort(c("sd_AVAL"), dirs = "desc") %>% 
   set_tags(
      stats = c("mean_AVAL","sd_AVAL","max_AVAL","mean_CHG","sd_CHG","max_CHG"),
      info  = c("USUBJID", "TRT01P")) %>% 
   view_trelliscope()

knitr::include_app('_trelliscope/index.html', height = "900px")

Scatterplot Matrix

Code 
pacman::p_load(scatterPlotMatrix)

d <- dat_twisk_w[,-c(1:2)] %>%
   relocate(TRT01P, .after = last_col()) %>% 
   na.omit()

scatterPlotMatrix(data = d,
                  zAxisDim = 'TRT01P',
                  categoricalCS = "Set1",
                  slidersPosition = list(
                     dimCount = 4,
                     xStartingDimIndex = 2,
                     yStartingDimIndex = 2 ),
                  columnLabels = c('W0', 'W12', 'W24','TRT'),
                  distribType = 1,
                  regressionType = 2,
                  plotProperties =  list( noCatColor= "#43665e", 
                                          watermarkColor = "#ddd", 
                                          point = list( alpha = 0.5, radius = 3 ),
                                          regression = list( strokeWidth = 1 ) ),
                  corrPlotType = "Text",
                  width = '1000px', height = '1000px')

Mixed models for repeated measures (MMRM)

Method 1: Longitudinal analysis of covariance

(1a) Overall TRT effect

\[Y_{t} = \beta_0 + \beta_1X + \beta_2Y_{t0}\]

Code 
m1a <- mmrm(AVAL ~ TRT01P + BASE + 
               us(AVISIT|USUBJID),     
            data = dat_twisk_l %>% filter(AVISITN != 0))
e1a <- emmeans(m1a, revpairwise ~ TRT01P) 
e1a
$emmeans
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       127 1.10 119      125      129 115.409  <.0001
 TRT       124 1.13 117      121      126 109.409  <.0001

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95 

$contrasts
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -3.69 1.59 118    -6.83   -0.548  -2.325  0.0218

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95
Code 
m1a$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
AVAL TRT01P BASE USUBJID AVISIT
1 129 TRT 124 S-00002 Week 12
2 126 TRT 124 S-00002 Week 24
3 132 CNT 119 S-00003 Week 12
4 114 CNT 119 S-00003 Week 24
5 115 CNT 113 S-00004 Week 12
6..248
249 118 TRT 131 S-00152 Week 24
Code 
tbl_regression(m1a,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 49 6.91 7.02 35, 62 <0.001
BASE 0.61 0.052 11.7 0.51, 0.72 <0.001
TRT01P




    CNT — — — —
    TRT -3.7 1.59 -2.33 -6.8, -0.55 0.022
AIC = 1,864; BIC = 1,873; Log-likelihood = -929; Deviance = 1,858
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
where AVISITN ^= 0;
class USUBJID AVISIT (ref = first) TRT01P (ref = first);
model AVAL = TRT01P BASE / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans TRT01P / diff cl;
run;

(1b) Weekly TRT effect

\[Y_{t} = \beta_0 + \beta_1X + \beta_2Y_{t0} + \beta_3time + \beta_4X\times time\]

Code 
m1b <- mmrm(AVAL ~ TRT01P + BASE + AVISIT + AVISIT:TRT01P + 
               us(AVISIT|USUBJID),     
            data = dat_twisk_l %>% filter(AVISITN != 0))
e1b <- emmeans(m1b, revpairwise ~ TRT01P | AVISIT) 
e1b
$emmeans
AVISIT = Week 12:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       129 1.30 127      126      131  98.565  <.0001
 TRT       124 1.34 127      121      127  92.235  <.0001

AVISIT = Week 24:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       126 1.35 115      123      128  93.294  <.0001
 TRT       123 1.36 114      120      126  90.192  <.0001

Confidence level used: 0.95 

$contrasts
AVISIT = Week 12:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -4.58 1.88 128    -8.30   -0.858  -2.434  0.0163

AVISIT = Week 24:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -2.68 1.92 115    -6.49    1.132  -1.393  0.1664

Confidence level used: 0.95
Code 
m1b$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
AVAL TRT01P BASE AVISIT USUBJID
1 129 TRT 124 Week 12 S-00002
2 126 TRT 124 Week 24 S-00002
3 132 CNT 119 Week 12 S-00003
4 114 CNT 119 Week 24 S-00003
5 115 CNT 113 Week 12 S-00004
6..248
249 118 TRT 131 Week 24 S-00152
Code 
tbl_regression(m1b,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 50 6.98 7.17 36, 64 <0.001
AVISIT




    Week 12 — — — —
    Week 12 — — — —
    Week 24 -2.9 1.46 -1.97 -5.8, 0.02 0.052
BASE 0.61 0.052 11.7 0.51, 0.72 <0.001
TRT01P




    CNT — — — —
    TRT -4.6 1.88 -2.43 -8.3, -0.86 0.016
TRT01P * AVISIT




    TRT * Week 24 1.9 2.08 0.914 -2.2, 6.0 0.36
AIC = 1,855; BIC = 1,863; Log-likelihood = -924; Deviance = 1,849
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
where AVISITN ^= 0;
class USUBJID AVISIT (ref = first) TRT01P (ref = first);
model AVAL = TRT01P BASE AVISIT AVISIT*TRT01P/ solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans AVISIT*TRT01P / slicediff=AVISIT diff cl;
run;

Method 2: Repeated measures

(2a) Overall TRT effect

\[Y_{t} = \beta_0 + \beta_1X + \beta_2time + \beta_3time\times X\]

Code 
m2a <- mmrm(AVAL ~ TRT01P + ABLFL + ABLFL:TRT01P +
               us(AVISIT|USUBJID),     
            data = dat_twisk_l)
e2a <- emmeans(m2a, revpairwise ~ TRT01P | ABLFL) 
e2a
$emmeans
ABLFL = N:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       128 1.55 133      125      132  82.619  <.0001
 TRT       122 1.59 132      119      126  77.020  <.0001

ABLFL = Y:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       130 1.80 137      127      134  72.239  <.0001
 TRT       126 1.84 137      123      130  68.592  <.0001

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95 

$contrasts
ABLFL = N:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -6.07 2.22 133   -10.47    -1.68  -2.732  0.0072

ABLFL = Y:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -3.89 2.58 137    -8.99     1.21  -1.508  0.1338

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95
Code 
m2a$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
AVAL TRT01P ABLFL USUBJID AVISIT
1 124 TRT Y S-00002 Week 0
2 129 TRT N S-00002 Week 12
3 126 TRT N S-00002 Week 24
4 119 CNT Y S-00003 Week 0
5 132 CNT N S-00003 Week 12
6..387
388 118 TRT N S-00152 Week 24
Code 
tbl_regression(m2a,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 128 1.55 82.6 125, 132 <0.001
ABLFL




    N — — — —
    Y 1.8 1.30 1.42 -0.73, 4.4 0.16
TRT01P




    CNT — — — —
    TRT -6.1 2.22 -2.73 -10, -1.7 0.007
TRT01P * ABLFL




    TRT * Y 2.2 1.86 1.17 -1.5, 5.9 0.24
AIC = 3,009; BIC = 3,027; Log-likelihood = -1,498; Deviance = 2,997
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
class USUBJID AVISIT (ref = first) ABLFL (ref = first) TRT01P (ref = first);
model AVAL =  TRT01P ABLFL ABLFL*TRT01P  / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans TRT01P*ABLFL / slicediff=ABLFL diff cl;
run;

(2b) Weekly TRT effect

\[Y_{t} = \beta_0 + \beta_1X + \beta_2dummy\_time_1 + \beta_3dummy\_time_2 +\\ \beta_4dumm\_time_1\times X + \beta_5dummy\_time_2 \times X\]

Code 
m2b <- mmrm(AVAL ~ TRT01P + AVISIT  + AVISIT :TRT01P +
               us(AVISIT|USUBJID),     
            data = dat_twisk_l)
e2b <- emmeans(m2b, revpairwise ~ TRT01P | AVISIT ) 
e2b
$emmeans
AVISIT = Week 0:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       131 1.81 137      127      134  72.003  <.0001
 TRT       127 1.85 137      123      130  68.229  <.0001

AVISIT = Week 12:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       130 1.76 134      127      134  74.083  <.0001
 TRT       123 1.80 135      119      127  68.204  <.0001

AVISIT = Week 24:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT       127 1.68 129      124      130  75.574  <.0001
 TRT       122 1.71 127      119      125  71.349  <.0001

Confidence level used: 0.95 

$contrasts
AVISIT = Week 0:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -4.15 2.59 137    -9.28    0.984  -1.598  0.1123

AVISIT = Week 12:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -7.09 2.52 134   -12.07   -2.117  -2.819  0.0055

AVISIT = Week 24:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -5.23 2.40 128    -9.97   -0.482  -2.180  0.0311

Confidence level used: 0.95
Code 
m2b$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
AVAL TRT01P AVISIT USUBJID
1 124 TRT Week 0 S-00002
2 129 TRT Week 12 S-00002
3 126 TRT Week 24 S-00002
4 119 CNT Week 0 S-00003
5 132 CNT Week 12 S-00003
6..387
388 118 TRT Week 24 S-00152
Code 
tbl_regression(m2b,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 131 1.81 72.0 127, 134 <0.001
AVISIT




    Week 0 — — — —
    Week 12 -0.63 1.44 -0.441 -3.5, 2.2 0.66
    Week 24 -3.5 1.56 -2.27 -6.6, -0.45 0.025
TRT01P




    CNT — — — —
    TRT -4.1 2.59 -1.60 -9.3, 0.98 0.11
TRT01P * AVISIT




    TRT * Week 12 -2.9 2.06 -1.43 -7.0, 1.1 0.16
    TRT * Week 24 -1.1 2.22 -0.486 -5.5, 3.3 0.63
AIC = 2,999; BIC = 3,017; Log-likelihood = -1,494; Deviance = 2,987
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
class USUBJID AVISIT (ref = first) TRT01P (ref = first);
model AVAL =  TRT01P AVISIT AVISIT*TRT01P  / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans TRT01P*AVISIT / slicediff=AVISIT diff cl;
run;

Method 2: Repeated measures wo/ treatment

(2c) Overall TRT effect

\[Y_{t} = \beta_0 + \beta_1time + \beta_2time\times X\]

Code 
m2c <- mmrm(AVAL ~ ABLFL + ABLFL:TRT01P +
               us(AVISIT|USUBJID),     
            data = dat_twisk_l %>% mutate(ABLFL = as.numeric(ABLFL)-2) )
e2cf0 <- emmeans(m2c,          ~ TRT01P | ABLFL, at = list( ABLFL = 0) ) 
e2cf1 <- emmeans(m2c, pairwise ~ TRT01P | ABLFL, at = list( ABLFL = 1) ) 
e2cf1
$emmeans
ABLFL = 1:
 TRT01P emmean SE df asymp.LCL asymp.UCL z.ratio p.value
 CNT    nonEst NA NA        NA        NA      NA      NA
 TRT    nonEst NA NA        NA        NA      NA      NA

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95 

$contrasts
ABLFL = 1:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 CNT - TRT     -3.7 1.57 118    -6.81    -0.59  -2.356  0.0201

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95
Code 
m2c$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
AVAL ABLFL USUBJID AVISIT TRT01P
1 124 0 S-00002 Week 0 TRT
2 129 -1 S-00002 Week 12 TRT
3 126 -1 S-00002 Week 24 TRT
4 119 0 S-00003 Week 0 CNT
5 132 -1 S-00003 Week 12 CNT
6..387
388 118 -1 S-00152 Week 24 TRT
Code 
tbl_regression(m2c,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 128 1.29 99.1 126, 131 <0.001
ABLFL 1.1 1.21 0.905 -1.3, 3.5 0.37
ABLFL * TRT01P




    ABLFL * TRT 3.7 1.57 2.36 0.59, 6.8 0.020
AIC = 3,015; BIC = 3,033; Log-likelihood = -1,501; Deviance = 3,003
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
class USUBJID TRT01P (ref = first);
model AVAL =  ABLFLN ABLFLN*TRT01P  / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
run;

(2d) Weekly TRT effect

\[Y_{t} = \beta_0 + \beta_1dummytime_1 + \beta_2dummytime_2 + \\ \beta_3dummytime_1 \times X + \beta_4dummytime_2 \times X\]

Code 
m2d <- mmrm(AVAL ~ I(AVISIT == "Week 12") + I(AVISIT == "Week 24") + 
               I(AVISIT == "Week 12" & TRT01P ==  "TRT") + I(AVISIT == "Week 24" & TRT01P == "TRT") + 
               us(AVISIT|USUBJID),     
            data = dat_twisk_l )

summary(m2d)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')
Estimate Std. Error t value Pr(>|t|)
(Intercept) 128.65 1.30 98.63 <0.001
I(AVISIT == "Week 12")TRUE   0.07 1.37  0.05  0.96 
I(AVISIT == "Week 24")TRUE  −2.64 1.46 −1.81  0.072
I(AVISIT == "Week 12" & TRT01P == "TRT")TRUE  −4.38 1.86 −2.36  0.020
I(AVISIT == "Week 24" & TRT01P == "TRT")TRUE  −2.90 1.91 −1.52  0.13 
Code 
m2d$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
AVAL I(AVISIT == "Week 12") I(AVISIT == "Week 24") I(AVISIT == "Week 12" & TRT01P == "TRT") I(AVISIT == "Week 24" & TRT01P == "TRT") USUBJID AVISIT
1 124 FALSE FALSE FALSE FALSE S-00002 Week 0
2 129 TRUE FALSE TRUE FALSE S-00002 Week 12
3 126 FALSE TRUE FALSE TRUE S-00002 Week 24
4 119 FALSE FALSE FALSE FALSE S-00003 Week 0
5 132 TRUE FALSE FALSE FALSE S-00003 Week 12
6..387
388 118 FALSE TRUE FALSE TRUE S-00152 Week 24
Code 
tbl_regression(m2d,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 129 1.30 98.6 126, 131 <0.001
I(AVISIT == "Week 12" & TRT01P == "TRT")




    FALSE — — — —
    TRUE -4.4 1.86 -2.36 -8.1, -0.71 0.020
I(AVISIT == "Week 12")




    FALSE — — — —
    TRUE 0.07 1.37 0.049 -2.6, 2.8 0.96
I(AVISIT == "Week 24" & TRT01P == "TRT")




    FALSE — — — —
    TRUE -2.9 1.91 -1.52 -6.7, 0.87 0.13
I(AVISIT == "Week 24")




    FALSE — — — —
    TRUE -2.6 1.46 -1.81 -5.5, 0.24 0.072
AIC = 3,006; BIC = 3,023; Log-likelihood = -1,497; Deviance = 2,994
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
class USUBJID TRT01P (ref = first);
model AVAL = Week_12 Week_24 Week_12*TRT01P Week_24*TRT01P / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
run;

Method 3: Analysis of changes (not adjusted)

(3a) Overall TRT effect

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X\]

Code 
m3a <- mmrm(CHG ~ TRT01P + 
               us(AVISIT|USUBJID),     
            data = dat_twisk_l %>% filter(AVISITN != 0))
e3a <- emmeans(m3a, revpairwise ~ TRT01P) 
e3a
$emmeans
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT     -1.85 1.31 120    -4.44     0.74  -1.415  0.1598
 TRT     -3.81 1.34 118    -6.47    -1.15  -2.835  0.0054

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95 

$contrasts
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -1.96 1.88 119    -5.67     1.75  -1.044  0.2984

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95
Code 
m3a$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
CHG TRT01P USUBJID AVISIT
1 5 TRT S-00002 Week 12
2 2 TRT S-00002 Week 24
3 13 CNT S-00003 Week 12
4 -5 CNT S-00003 Week 24
5 2 CNT S-00004 Week 12
6..248
249 -13 TRT S-00152 Week 24
Code 
tbl_regression(m3a,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) -1.9 1.31 -1.41 -4.4, 0.74 0.16
TRT01P




    CNT — — — —
    TRT -2.0 1.88 -1.04 -5.7, 1.8 0.30
AIC = 1,906; BIC = 1,914; Log-likelihood = -950; Deviance = 1,900
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
where AVISITN ^= 0;
class USUBJID TRT01P (ref = first);
model CHG = TRT01P / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans TRT01P / diff cl;
run;

(3b) Weekly TRT effect

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X + \beta_2 time + \beta_3X \times time\]

Code 
m3b <- mmrm(CHG ~ TRT01P + AVISIT + TRT01P:AVISIT +
               us(AVISIT|USUBJID),     
            data = dat_twisk_l %>% filter(AVISITN != 0))
e3b <- emmeans(m3b, revpairwise ~ TRT01P | AVISIT) 
e3b
$emmeans
AVISIT = Week 12:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT    -0.642 1.44 128    -3.50    2.216  -0.444  0.6575
 TRT    -3.381 1.49 128    -6.33   -0.434  -2.270  0.0249

AVISIT = Week 24:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT    -3.563 1.57 112    -6.68   -0.447  -2.265  0.0254
 TRT    -4.379 1.60 110    -7.54   -1.216  -2.744  0.0071

Confidence level used: 0.95 

$contrasts
AVISIT = Week 12:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT   -2.739 2.07 128    -6.84     1.37  -1.320  0.1891

AVISIT = Week 24:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT   -0.817 2.24 111    -5.26     3.62  -0.364  0.7162

Confidence level used: 0.95
Code 
m3a$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
CHG TRT01P USUBJID AVISIT
1 5 TRT S-00002 Week 12
2 2 TRT S-00002 Week 24
3 13 CNT S-00003 Week 12
4 -5 CNT S-00003 Week 24
5 2 CNT S-00004 Week 12
6..248
249 -13 TRT S-00152 Week 24
Code 
tbl_regression(m3b,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) -0.64 1.44 -0.444 -3.5, 2.2 0.66
AVISIT




    Week 12 — — — —
    Week 12 — — — —
    Week 24 -2.9 1.46 -1.99 -5.8, -0.02 0.048
TRT01P




    CNT — — — —
    TRT -2.7 2.07 -1.32 -6.8, 1.4 0.19
TRT01P * AVISIT




    TRT * Week 24 1.9 2.08 0.923 -2.2, 6.1 0.36
AIC = 1,896; BIC = 1,905; Log-likelihood = -945; Deviance = 1,890
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
where AVISITN ^= 0;
class USUBJID AVISIT (ref = first) TRT01P (ref = first);
model CHG = TRT01P AVISIT AVISIT*TRT01P / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans TRT01P*AVISIT / slicediff=AVISIT diff cl;
run;

Method 3: Analysis of changes (adjusted)

(3c) Overall TRT effect

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X +\beta_2Y_{t0}\]

Code 
m3c <- mmrm(CHG ~ TRT01P + BASE +
               us(AVISIT|USUBJID),     
            data = dat_twisk_l %>% filter(AVISITN != 0))
e3c <- emmeans(m3c, revpairwise ~ TRT01P) 
e3c
$emmeans
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT     -1.06 1.10 119    -3.24     1.12  -0.963  0.3373
 TRT     -4.75 1.13 117    -6.99    -2.52  -4.210  0.0001

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95 

$contrasts
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -3.69 1.59 118    -6.83   -0.548  -2.325  0.0218

Results are averaged over the levels of: AVISIT 
Confidence level used: 0.95
Code 
m3c$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
CHG TRT01P BASE USUBJID AVISIT
1 5 TRT 124 S-00002 Week 12
2 2 TRT 124 S-00002 Week 24
3 13 CNT 119 S-00003 Week 12
4 -5 CNT 119 S-00003 Week 24
5 2 CNT 113 S-00004 Week 12
6..248
249 -13 TRT 131 S-00152 Week 24
Code 
tbl_regression(m3c,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic)) %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 49 6.91 7.02 35, 62 <0.001
BASE -0.39 0.052 -7.40 -0.49, -0.28 <0.001
TRT01P




    CNT — — — —
    TRT -3.7 1.59 -2.33 -6.8, -0.55 0.022
AIC = 1,864; BIC = 1,873; Log-likelihood = -929; Deviance = 1,858
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
where AVISITN ^= 0;
class USUBJID TRT01P (ref = first);
model CHG = TRT01P BASE / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans TRT01P / diff cl;
run;

(3d) Weekly TRT effect

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X + \beta_2 Y_{t0} + \beta_3X \times time\]

Code 
m3d <- mmrm(CHG ~ TRT01P + BASE + AVISIT + TRT01P:AVISIT +
               us(AVISIT|USUBJID),     
            data = dat_twisk_l %>% filter(AVISITN != 0))
e3d <- emmeans(m3d, revpairwise ~ TRT01P | AVISIT) 
e3d
$emmeans
AVISIT = Week 12:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT      0.30 1.30 127    -2.28   2.8808   0.230  0.8187
 TRT     -4.28 1.34 127    -6.94  -1.6206  -3.184  0.0018

AVISIT = Week 24:
 TRT01P emmean   SE  df lower.CL upper.CL t.ratio p.value
 CNT     -2.57 1.35 115    -5.24   0.0969  -1.909  0.0588
 TRT     -5.25 1.36 114    -7.95  -2.5502  -3.851  0.0002

Confidence level used: 0.95 

$contrasts
AVISIT = Week 12:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -4.58 1.88 128    -8.30   -0.858  -2.434  0.0163

AVISIT = Week 24:
 contrast  estimate   SE  df lower.CL upper.CL t.ratio p.value
 TRT - CNT    -2.68 1.92 115    -6.49    1.132  -1.393  0.1664

Confidence level used: 0.95
Code 
m3d$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)
CHG TRT01P BASE AVISIT USUBJID
1 5 TRT 124 Week 12 S-00002
2 2 TRT 124 Week 24 S-00002
3 13 CNT 119 Week 12 S-00003
4 -5 CNT 119 Week 24 S-00003
5 2 CNT 113 Week 12 S-00004
6..248
249 -13 TRT 131 Week 24 S-00152
Code 
tbl_regression(m3d,
               intercept = TRUE,
               pvalue_fun = function(x) style_pvalue(x, digits = 2) ) %>% 
   modify_column_unhide(column = c(std.error, statistic))  %>% 
   add_glance_source_note()
Characteristic Beta SE1 Statistic 95% CI1 p-value
(Intercept) 50 6.98 7.17 36, 64 <0.001
AVISIT




    Week 12 — — — —
    Week 12 — — — —
    Week 24 -2.9 1.46 -1.97 -5.8, 0.02 0.052
BASE -0.39 0.052 -7.39 -0.49, -0.28 <0.001
TRT01P




    CNT — — — —
    TRT -4.6 1.88 -2.43 -8.3, -0.86 0.016
TRT01P * AVISIT




    TRT * Week 24 1.9 2.08 0.914 -2.2, 6.0 0.36
AIC = 1,855; BIC = 1,863; Log-likelihood = -924; Deviance = 1,849
1 SE = Standard Error, CI = Confidence Interval
Code 
proc glimmix data=dat_twisk_l noclprint = 10;
where AVISITN ^= 0;
class USUBJID AVISIT (ref = first) TRT01P (ref = first);
model CHG = TRT01P BASE AVISIT AVISIT*TRT01P / solution ddfm=satterthwaite;
random _residual_ / subject=USUBJID type=un;
lsmeans TRT01P*AVISIT / slicediff=AVISIT diff cl;
run;

Overall TRT Effect (Combined Results)

Beta SE Statistic 95% CI p-value model
(1a) Longitudinal analysis of covariance −3.7 1.59 −2.33  (−6.8, −0.5)  0.022 AVAL ~ TRT01P + BASE + us(AVISIT | USUBJID)
(2a) Repeated measures analysis −6.1 2.22 −2.73 (−10.5, −1.7)  0.007 AVAL ~ TRT01P + ABLFL + ABLFL:TRT01P + us(AVISIT | USUBJID)
(2c) Repeated measures wo/ treatment −3.7 1.57 −2.36  (−6.8, −0.6)  0.020 AVAL ~ ABLFL + ABLFL:TRT01P + us(AVISIT | USUBJID)
(3a) Analysis of changes (not adjusted) −2.0 1.88 −1.04  (−5.7, 1.8)   0.30  CHG ~ TRT01P + us(AVISIT | USUBJID)
(3c) Analysis of changes (adjusted) −3.7 1.59 −2.33  (−6.8, −0.5)  0.022 CHG ~ TRT01P + BASE + us(AVISIT | USUBJID)
  • (1a) overall TRT effect of -3.7 (mmHg)
  • (2a) wo/ adjustment for baseline differences overestimation
  • (2b) w/ adjustment for baseline differences similar results to (1a)
  • (3a) wo/ adjustments for baseline difference underestimation
  • (3c) w/ adjustment for baseline differences exact same results

Mathematical equivalence between models


(1a) Longitudinal analysis of covariance (3c) Analysis of changes (adjusted)
\[Y_{t} = \beta_0 + \beta_1X + \beta_2Y_{t0}\] \[Y_{t} - Y_{t0} = \beta_0 + \beta_1X +\beta_2Y_{t0}\]
\[Y_{t} = \beta_0 + \beta_1X +\beta_2Y_{t0} + Y_{t0}\]
\[Y_{t} = \beta_0 + \beta_1X + (1 + \beta_2) Y_{t0}\]
\[Y_{t} = 49 - 3.7X + 0.6Y_{t0}\] \[Y_{t} = 49 - 3.7X - 0.4Y_{t0}\]

ADaM Basic Data Structure (BDS)

Simulation (Complete Data)

Code 
sim_1 <- import('mod/sim/sim_1.rds')

sim_w <- sim_1 %>% 
   select(-BASE, -CHG) %>% 
   pivot_wider(names_from = c(MISS, AVISIT),
               values_from = c(AVAL) )

add_by_n <- function(data, variable, by, tbl, ...) {
   data %>%
      select(all_of(c(variable, by))) %>%
      dplyr::group_by(.data[[by]]) %>%
      dplyr::summarise_all(~sum(!is.na(.))) %>%
      rlang::set_names(c("by", "variable")) %>%
      dplyr::left_join(
         tbl$df_by %>% select(by, by_col),
         by = "by"
      ) %>%
      mutate(
         by_col = paste0("add_n_", by_col),
         variable = style_number(variable)
      ) %>%
      select(-by) %>%
      tidyr::pivot_wider(names_from = by_col, 
                         values_from = variable)
}

tbl_summary(data = sim_w %>% mutate(TRT01P = fct_rev(TRT01P)),
            by = TRT01P , 
            include = starts_with('COMP_'),
            digit = list(everything() ~ c(1,1)),
            statistic = list(all_continuous() ~ "{mean} ({sd})"),
            missing = 'no',
            label = list( `COMP_Week 0` ~ 'Week 0 (Randomization)',
                          `COMP_Week 12` ~ 'Week 12',
                          `COMP_Week 24` ~ 'Week 24')) %>% 
   add_stat( fns = everything() ~ add_by_n ) %>%
   add_difference(test = list(all_continuous() ~ "t.test"),
                  pvalue_fun = ~style_pvalue(.x, digits = 2),
                  estimate_fun = list(all_continuous() ~ function(x) style_number(x, digits=1)) ) %>% 
   bold_labels() %>% 
   modify_column_merge(pattern = "{estimate} ({ci})") %>% 
   modify_header(
      label = '',
      estimate = '**DIFF (95% CI)**',
      all_stat_cols(FALSE) ~ "**{level}**",
      starts_with("add_n_stat") ~ "**n**") %>% 
   modify_table_body(
      ~ .x %>%
         dplyr::relocate(add_n_stat_1, .before = stat_1) %>%
         dplyr::relocate(add_n_stat_2, .before = stat_2)
   ) %>% 
   modify_spanning_header(
      c(all_stat_cols(F), starts_with("add_n_stat")) ~ "**Treatment**") %>% 
   as_gt() %>% 
   gtsave('tbl/sim-tbl_comp.png', zoom = 3)

knitr::include_graphics('tbl/sim-tbl_comp.png')

Simulation (Missing Data)

Code 
tbl_summary(data = sim_w %>% mutate(TRT01P = fct_rev(TRT01P)),
            by = TRT01P , 
            include = starts_with('MISS_'),
            digit = list(everything() ~ c(1,1)),
            statistic = list(all_continuous() ~ "{mean} ({sd})"),
            missing = 'no',
            label = list( `MISS_Week 0` ~ 'Week 0 (Randomization)',
                          `MISS_Week 12` ~ 'Week 12',
                          `MISS_Week 24` ~ 'Week 24')) %>% 
   add_stat( fns = everything() ~ add_by_n ) %>%
   add_difference(test = list(all_continuous() ~ "t.test"),
                  pvalue_fun = ~style_pvalue(.x, digits = 2),
                  estimate_fun = list(all_continuous() ~ function(x) style_number(x, digits=1)) ) %>% 
   bold_labels() %>% 
   modify_column_merge(pattern = "{estimate} ({ci})") %>% 
   modify_header(
      label = '',
      estimate = '**DIFF (95% CI)**',
      all_stat_cols(FALSE) ~ "**{level}**",
      starts_with("add_n_stat") ~ "**n**") %>% 
   modify_table_body(
      ~ .x %>%
         dplyr::relocate(add_n_stat_1, .before = stat_1) %>%
         dplyr::relocate(add_n_stat_2, .before = stat_2)
   ) %>% 
   modify_spanning_header(
      c(all_stat_cols(F), starts_with("add_n_stat")) ~ "**Treatment**") %>% 
   as_gt() %>% 
   gtsave('tbl/sim-tbl_miss.png', zoom = 3)

knitr::include_graphics('tbl/sim-tbl_miss.png')

Simulation (Power Analyses)

Code 
pacman::p_load(ggdist)

sim_l <- import('mod/sim/sim_l.rds')

sim_p <- sim_l %>% 
   group_by(model, MISS) %>% 
   summarise(power = sum(I(statistic< -1.98)) / n() )

sim_g <- sim_l %>% 
   group_by(model, MISS) %>% 
   summarise(statistic = mean(statistic) )

sim_l %>%
   ggplot(aes(y = model, x = statistic )) +
   geom_vline(xintercept  = -2, color = 'gray50') +
   stat_halfeye(adjust = 0.70, alpha = 0.50) +
   facet_wrap(~ MISS,
              labeller = as_labeller(c(COMP='Completed Data', MISS='MCAR (25%)')) ) +
   scale_y_discrete(limits = rev,
                    name = '') +
   geom_text(data = sim_p,
             vjust = -1.0,
             aes(x = 0.85, label = style_percent(power, digits = 0, symbol = TRUE) ) ) +
   geom_text(data = sim_g,
             vjust = -1.0,
             aes(x = statistic, label = style_number(statistic, digits = 1) ) ) +
   theme_bw(base_size = 18) +
   theme(panel.grid.minor = element_blank(),
         axis.text.y = element_text(face = 'bold',size = 10))

Different ways to estimate treatment effects in randomised controlled trials

Data Visualizations for Clinical Trials Reporting
Some Examples
Visualizations for the Special Interest Group at
Statisticians in the Pharmaceutical Industry (PSI)

Quality of life outcomes in a cancer trial

Patient Profiler (Proof of Concept)

Sensitivity Analyses of HiSCR definition on the results

Relative Importance of Regressors in Linear Models

Continuous Glucose Monitoring (CGM) Visualization

Prediction of Residual Tumor model w/ ModelStudio

Competing Risk Tables Validation Visualization