MMRM table and plots

Author

Affiliation

Eric

test

Code

rm(list=ls(all=TRUE))
library(tinytex)
library(ggplot2)
library(plotly)
library(dplyr)
library(haven)
library(labelled)
library(bslib)
library(DT)
library(r2rtf)
library(metalite)
library(metalite.table1)
library(haven)
library(tidyr)
library(tools)
library(gtsummary)
library(stringr)
library(RColorBrewer)
library(flextable)
library(stringr)
library(tern)
library(tidyverse)
library(vistime)
library(readxl)
library(hrbrthemes)
library(rlang)
library(lubridate) # For date calculations
library(haven)
library(metalite)
library(metalite.table1)

library(tern.mmrm)
library(broom)
library(random.cdisc.data)
library(MASS)

Code

adsl <- random.cdisc.data::cadsl
adqs <- random.cdisc.data::cadqs

# Ensure character variables are converted to factors and empty strings and NAs are explicit missing levels.
adsl <- df_explicit_na(adsl)
adqs <- df_explicit_na(adqs)

adqs_f <- adqs %>%
  dplyr::filter(PARAMCD == "FKSI-FWB" & !AVISIT %in% c("BASELINE", "SCREENING")) %>%
  droplevels() %>%
  dplyr::mutate(ARMCD = factor(ARMCD, levels = c("ARM B", "ARM A", "ARM C"))) %>%
  dplyr::mutate(
    AVISITN = rank(AVISITN) %>%
      as.factor() %>%
      as.numeric() %>%
      as.factor()
  )
adsl_sub <- adqs_f %>%
  dplyr::filter(!is.na(CHG)) %>%
  distinct(USUBJID) %>%
  left_join(adsl, by = "USUBJID")
var_labels(adqs_f) <- var_labels(adqs)

Code

MASS::truehist(adqs_f$CHG, nbins = 12)

Code

mmrm_results <- fit_mmrm(
  vars = list(
    response = "CHG",
    covariates = c("BASE", "STRATA1", "BMRKR2"),
    id = "USUBJID",
    arm = "ARMCD",
    visit = "AVISIT"
  ),
  data = adqs_f,
  weights_emmeans = "equal"
)

mmrm() registered as emmeans extension

Code

df <- tidy(mmrm_results)
attr(df$AVISIT, "label") <- "Visit"

# Define the split function
split_fun <- drop_split_levels

result <- basic_table(show_colcounts = TRUE) %>%
  split_cols_by("ARMCD", ref_group = mmrm_results$ref_level) %>%
  split_rows_by("AVISIT", split_fun = split_fun, label_pos = "topleft", split_label = obj_label(df$AVISIT)) %>%
  summarize_lsmeans(show_relative = "increase") %>%
  append_topleft("  Statistics") %>%
  build_table(df, alt_counts_df = adsl_sub)

result

Visit                                      ARM B              ARM A              ARM C      
  Statistics                              (N=134)            (N=134)            (N=132)     
--------------------------------------------------------------------------------------------
WEEK 1 DAY 8                                                                                
  n                                         134                134                132       
  Adjusted Mean (SE)                   3.488 (0.687)      4.246 (0.687)      3.163 (0.692)  
    95% CI                             (2.136, 4.839)     (2.895, 5.598)     (1.803, 4.523) 
  Difference in Adjusted Means (SE)                       0.759 (0.973)      -0.325 (0.976) 
    95% CI                                               (-1.155, 2.672)    (-2.243, 1.593) 
    Relative Increase (%)                                     21.8%              -9.3%      
    p-value (MMRM)                                            0.4362             0.7394     
WEEK 2 DAY 15                                                                               
  n                                         134                134                132       
  Adjusted Mean (SE)                   9.135 (0.768)      9.018 (0.767)      8.509 (0.773)  
    95% CI                            (7.626, 10.644)    (7.510, 10.527)    (6.991, 10.028) 
  Difference in Adjusted Means (SE)                       -0.117 (1.087)     -0.626 (1.090) 
    95% CI                                               (-2.253, 2.020)    (-2.768, 1.517) 
    Relative Increase (%)                                     -1.3%              -6.8%      
    p-value (MMRM)                                            0.9146             0.5662     
WEEK 3 DAY 22                                                                               
  n                                         134                134                132       
  Adjusted Mean (SE)                   13.547 (0.871)     16.014 (0.871)     15.789 (0.877) 
    95% CI                            (11.835, 15.259)   (14.303, 17.726)   (14.066, 17.513)
  Difference in Adjusted Means (SE)                       2.467 (1.232)      2.242 (1.236)  
    95% CI                                                (0.044, 4.890)    (-0.187, 4.672) 
    Relative Increase (%)                                     18.2%              16.6%      
    p-value (MMRM)                                            0.0460             0.0704     
WEEK 4 DAY 29                                                                               
  n                                         134                134                132       
  Adjusted Mean (SE)                   18.102 (0.995)     19.479 (0.995)     19.511 (1.002) 
    95% CI                            (16.146, 20.059)   (17.523, 21.435)   (17.541, 21.481)
  Difference in Adjusted Means (SE)                       1.377 (1.408)      1.409 (1.413)  
    95% CI                                               (-1.392, 4.145)    (-1.368, 4.186) 
    Relative Increase (%)                                      7.6%               7.8%      
    p-value (MMRM)                                            0.3288             0.3192     
WEEK 5 DAY 36                                                                               
  n                                         134                134                132       
  Adjusted Mean (SE)                   23.503 (1.058)     24.931 (1.058)     23.704 (1.065) 
    95% CI                            (21.423, 25.583)   (22.851, 27.011)   (21.610, 25.798)
  Difference in Adjusted Means (SE)                       1.428 (1.497)      0.201 (1.502)  
    95% CI                                               (-1.515, 4.371)    (-2.752, 3.153) 
    Relative Increase (%)                                      6.1%               0.9%      
    p-value (MMRM)                                            0.3408             0.8937     

Code

plot <- g_mmrm_lsmeans(
  mmrm_results,
  select = "estimates",
  xlab = "Visit"
)
plot

Code

plot <- g_mmrm_diagnostic(mmrm_results, type = "q-q-residual")
plot

Code

library(lme4)

Warning: package 'lme4' was built under R version 4.2.3

Loading required package: Matrix

Warning: package 'Matrix' was built under R version 4.2.3

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

Code

model <- lmer(CHG ~ AVISIT + ARMCD + (1|USUBJID), data = adqs_f)
library(broom.mixed)

Warning: package 'broom.mixed' was built under R version 4.2.3

Code

model_summary <- tidy(model, effects = "fixed")
model_summary

# A tibble: 7 x 5
  effect term                estimate std.error statistic
  <chr>  <chr>                  <dbl>     <dbl>     <dbl>
1 fixed  (Intercept)            3.64      0.972     3.75 
2 fixed  AVISITWEEK 2 DAY 15    5.25      0.729     7.21 
3 fixed  AVISITWEEK 3 DAY 22   11.5       0.729    15.8  
4 fixed  AVISITWEEK 4 DAY 29   15.4       0.729    21.1  
5 fixed  AVISITWEEK 5 DAY 36   20.4       0.729    28.0  
6 fixed  ARMCDARM A            -0.490     1.21     -0.405
7 fixed  ARMCDARM C             0.502     1.21      0.413

Code

library(emmeans)

Warning: package 'emmeans' was built under R version 4.2.3

Code

emmeans_results <- emmeans(model, specs = pairwise ~ ARMCD | AVISIT)
emmeans_summary <- summary(emmeans_results)
emmeans_summary

$emmeans
AVISIT = WEEK 1 DAY 8:
 ARMCD emmean    SE  df lower.CL upper.CL
 ARM B   3.64 0.972 647     1.73     5.55
 ARM A   3.15 0.972 647     1.24     5.06
 ARM C   4.14 0.978 643     2.22     6.06

AVISIT = WEEK 2 DAY 15:
 ARMCD emmean    SE  df lower.CL upper.CL
 ARM B   8.90 0.972 647     6.99    10.80
 ARM A   8.41 0.972 647     6.50    10.31
 ARM C   9.40 0.978 643     7.48    11.32

AVISIT = WEEK 3 DAY 22:
 ARMCD emmean    SE  df lower.CL upper.CL
 ARM B  15.12 0.972 647    13.21    17.03
 ARM A  14.63 0.972 647    12.72    16.54
 ARM C  15.62 0.978 643    13.70    17.54

AVISIT = WEEK 4 DAY 29:
 ARMCD emmean    SE  df lower.CL upper.CL
 ARM B  19.03 0.972 647    17.13    20.94
 ARM A  18.54 0.972 647    16.64    20.45
 ARM C  19.54 0.978 643    17.62    21.46

AVISIT = WEEK 5 DAY 36:
 ARMCD emmean    SE  df lower.CL upper.CL
 ARM B  24.05 0.972 647    22.14    25.96
 ARM A  23.56 0.972 647    21.65    25.47
 ARM C  24.56 0.978 643    22.64    26.47

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
AVISIT = WEEK 1 DAY 8:
 contrast      estimate   SE  df t.ratio p.value
 ARM B - ARM A    0.490 1.21 397   0.405  0.9137
 ARM B - ARM C   -0.502 1.21 397  -0.413  0.9103
 ARM A - ARM C   -0.991 1.21 397  -0.816  0.6933

AVISIT = WEEK 2 DAY 15:
 contrast      estimate   SE  df t.ratio p.value
 ARM B - ARM A    0.490 1.21 397   0.405  0.9137
 ARM B - ARM C   -0.502 1.21 397  -0.413  0.9103
 ARM A - ARM C   -0.991 1.21 397  -0.816  0.6933

AVISIT = WEEK 3 DAY 22:
 contrast      estimate   SE  df t.ratio p.value
 ARM B - ARM A    0.490 1.21 397   0.405  0.9137
 ARM B - ARM C   -0.502 1.21 397  -0.413  0.9103
 ARM A - ARM C   -0.991 1.21 397  -0.816  0.6933

AVISIT = WEEK 4 DAY 29:
 contrast      estimate   SE  df t.ratio p.value
 ARM B - ARM A    0.490 1.21 397   0.405  0.9137
 ARM B - ARM C   -0.502 1.21 397  -0.413  0.9103
 ARM A - ARM C   -0.991 1.21 397  -0.816  0.6933

AVISIT = WEEK 5 DAY 36:
 contrast      estimate   SE  df t.ratio p.value
 ARM B - ARM A    0.490 1.21 397   0.405  0.9137
 ARM B - ARM C   -0.502 1.21 397  -0.413  0.9103
 ARM A - ARM C   -0.991 1.21 397  -0.816  0.6933

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 3 estimates 

Code

 #Summarize the emmeans results for differences
emm_summary <- summary(pairs(emmeans_results), infer = c(TRUE, TRUE))

Code

library(gtsummary)
# Use gtsummary to create a table (gtsummary works well with data frames)
emm_df <- as.data.frame(emm_summary)

table <- emm_df %>%
  tbl_summary(
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ "{mean} ({sd})",
    digits = all_continuous() ~ 2
  ) %>%
  modify_header(label = "**Statistic**")

# Print the table with kableExtra for better formatting in R Markdown or output
library(kableExtra)

Warning: package 'kableExtra' was built under R version 4.2.3

Attaching package: 'kableExtra'

The following objects are masked from 'package:flextable':

    as_image, footnote

The following object is masked from 'package:dplyr':

    group_rows

Code

kable(table)

**Statistic**	**N = 15**
contrast	NA
ARM B - ARM A	5 (33%)
ARM B - ARM C	5 (33%)
ARM A - ARM C	5 (33%)
AVISIT	NA
WEEK 1 DAY 8	3 (20%)
WEEK 2 DAY 15	3 (20%)
WEEK 3 DAY 22	3 (20%)
WEEK 4 DAY 29	3 (20%)
WEEK 5 DAY 36	3 (20%)
estimate	NA
-0.991451914714403	5 (33%)
-0.501657129753606	5 (33%)
0.489794784960797	5 (33%)
SE	NA
1.21028619541399	5 (33%)
1.2148619629459	10 (67%)
df	NA
396.999999999982	10 (67%)
396.999999999994	5 (33%)
lower.CL	NA
-3.8494596648039	5 (33%)
-3.35966487984465	5 (33%)
-2.3574483025861	5 (33%)
upper.CL	NA
1.86655583537509	5 (33%)
2.35635062033744	5 (33%)
3.3370378725077	5 (33%)
t.ratio	NA
-0.816102524364373	5 (33%)
-0.412933440221591	5 (33%)
0.404693358328572	5 (33%)
p.value	NA
0.693263248485727	5 (33%)
0.910325443763207	5 (33%)
0.913709138581634	5 (33%)

Code

library(dplyr)
library(tern.mmrm)
library(nestcolor)

Warning: package 'nestcolor' was built under R version 4.2.3

Code

adsl <- random.cdisc.data::cadsl
adqs <- random.cdisc.data::cadqs

adqs_f <- adqs %>%
  dplyr::filter(PARAMCD == "FKSI-FWB" & !AVISIT %in% c("BASELINE")) %>%
  droplevels() %>%
  dplyr::mutate(ARMCD = factor(ARMCD, levels = c("ARM B", "ARM A", "ARM C"))) %>%
  dplyr::mutate(
    AVISITN = rank(AVISITN) %>%
      as.factor() %>%
      as.numeric() %>%
      as.factor()
  )

# Simulation of groups.
set.seed(2)
adqs_f_with_groups <- rbind(
  within(
    adqs_f[sample(seq_len(nrow(adqs_f)), size = 1 / 2 * nrow(adqs_f)), ],
    group <- "subpopulation 1"
  ),
  within(
    adqs_f,
    {
      group <- "subpopulation 2"
      AVAL <- AVAL + rnorm(length(AVAL), mean = 10, sd = 2)
      USUBJID <- paste0(USUBJID, "-S2")
    }
  )
)
adqs_f_with_groups$group <- factor(adqs_f_with_groups$group)


mmrm_results <- fit_mmrm(
  data = adqs_f_with_groups,
  vars = list(
    response = "AVAL",
    covariates = c(),
    id = "USUBJID",
    arm = "ARMCD",
    visit = "AVISIT"
  ),
  cor_struct = "unstructured",
  weights_emmeans = "equal",
  parallel = TRUE
)

df_a <- extract_mmrm_subgroups(
  fit = mmrm_results,
  visit = "SCREENING",
  subgroups = c("group", "SEX"),
  treatment_arm = "ARM A"
)

t<-tidy(mmrm_results)

df_a_t<-as.data.frame(df_a)
tab_a <- basic_table() %>%
  tabulate_mmrm_subgroups(
    df = df_a,
    vars = c("n_tot", "diff", "ci", "pval")
  )

plot <- g_forest(
  tab_a,
  logx = FALSE,
  xlim = c(-5, 2.5),
  x_at = c(-5, -2.5, 0, 2.5),
  vline = 0
)
plot

Code

path<-"//dc1/project workspace/DTX401 (GSD1a)/DTX401-CL301/11.0 Biostatistics/Programs/Dev/PgmStat/EC/SAS outputs/"

totaldose<-read_sas(paste0(path, "totaldose.sas7bdat"), NULL)
subj<-totaldose%>%filter(AVISITN==0)
subj$TRT01P<-factor(subj$TRT01P, levels = c("Placebo","DTX401"))

totaldose_f <- totaldose %>%
  dplyr::filter(AVISITN %in% c(168,210,252,294,336)) %>%

  dplyr::mutate(TRT01P = factor(TRT01P, levels = c("Placebo","DTX401"))) %>%
  dplyr::mutate(
    AVISITN = rank(AVISITN) %>%
      as.factor() %>%
      as.numeric() %>%
      as.factor()
  )

totaldose_f$AVISIT <- as.factor(totaldose_f$AVISIT)

Code

mmrm_results <- fit_mmrm(
  vars = list(
    response = "CHG",
    covariates = c("TRT01PN", "AVISITN" ,"TRT01PN:AVISITN", "STRATAGN", "STRATGLN",  "BASE"),
    id = "USUBJID",
     arm = "TRT01P",
    visit = "AVISIT"
  ),
  data = totaldose_f,
  weights_emmeans = "equal"
)

NOTE: A nesting structure was detected in the fitted model:
    TRT01PN %in% TRT01P, AVISITN %in% AVISIT

Code

df <- tidy(mmrm_results)
attr(df$AVISIT, "label") <- "Visit"

# Define the split function
split_fun <- drop_split_levels

result <- basic_table(show_colcounts = TRUE) %>%
  split_cols_by("TRT01P", ref_group = mmrm_results$ref_level) %>%
  split_rows_by("AVISIT", split_fun = split_fun, label_pos = "topleft", split_label = obj_label(df$AVISIT)) %>%
  summarize_lsmeans(show_relative = "increase") %>%
  append_topleft("  Statistics") %>%
  build_table(df, alt_counts_df = subj)

result

Visit                                     Placebo            DTX401     
  Statistics                              (N=24)             (N=20)     
------------------------------------------------------------------------
Week 24                                                                 
  n                                         24                 20       
  Adjusted Mean (SE)                  -0.236 (0.167)     -0.829 (0.185) 
    95% CI                            (-0.573, 0.101)   (-1.203, -0.455)
  Difference in Adjusted Means (SE)                      -0.593 (0.250) 
    95% CI                                              (-1.097, -0.088)
    Relative Increase (%)                                    251.1%     
    p-value (MMRM)                                           0.0226     
Week 30                                                                 
  n                                         23                 18       
  Adjusted Mean (SE)                  -0.018 (0.149)     -1.134 (0.166) 
    95% CI                            (-0.320, 0.283)   (-1.471, -0.798)
  Difference in Adjusted Means (SE)                      -1.116 (0.223) 
    95% CI                                              (-1.569, -0.663)
    Relative Increase (%)                                   6073.1%     
    p-value (MMRM)                                          <0.0001     
Week 36                                                                 
  n                                         22                 19       
  Adjusted Mean (SE)                   0.008 (0.141)     -1.014 (0.156) 
    95% CI                            (-0.278, 0.293)   (-1.330, -0.698)
  Difference in Adjusted Means (SE)                      -1.022 (0.211) 
    95% CI                                              (-1.449, -0.595)
    Relative Increase (%)                                  -13103.4%    
    p-value (MMRM)                                          <0.0001     
Week 42                                                                 
  n                                         22                 18       
  Adjusted Mean (SE)                  -0.214 (0.240)     -1.373 (0.266) 
    95% CI                            (-0.699, 0.271)   (-1.910, -0.836)
  Difference in Adjusted Means (SE)                      -1.159 (0.358) 
    95% CI                                              (-1.883, -0.435)
    Relative Increase (%)                                    541.7%     
    p-value (MMRM)                                           0.0024     
Week 48                                                                 
  n                                         23                 19       
  Adjusted Mean (SE)                  -0.169 (0.173)     -1.126 (0.192) 
    95% CI                            (-0.520, 0.183)   (-1.515, -0.737)
  Difference in Adjusted Means (SE)                      -0.957 (0.260) 
    95% CI                                              (-1.483, -0.432)
    Relative Increase (%)                                    568.1%     
    p-value (MMRM)                                           0.0007     

Code

plot <- g_mmrm_lsmeans(
  mmrm_results,
  select = "estimates",
  xlab = "Visit"
)
ggplotly(plot)

Code

plot <- g_mmrm_diagnostic(mmrm_results, type = "q-q-residual")
plot

Code

df_a <- extract_mmrm_subgroups(
  fit = mmrm_results,
  visit = "Week 48"
)

t<-tidy(mmrm_results)

df_a_t<-as.data.frame(df_a)
tab_a <- basic_table() %>%
  tabulate_mmrm_subgroups(
    df = df_a,
    vars = c("n_tot", "diff", "ci", "pval")
  )

plot <- g_forest(
  tab_a,
  logx = FALSE,
  xlim = c(-2, 1),
  x_at = c(-2, -1, 0, 1),
  vline = 0
)
plot

Citation

BibTeX citation:

@online{eric2024,
  author = {Eric},
  title = {MMRM Table and Plots},
  date = {2024-05-31},
  langid = {en}
}

For attribution, please cite this work as:

Eric. 2024. “MMRM Table and Plots.” Analysis. May 31, 2024.