code complete for emmeans contrasts helpers

Author: Daniel Sabanes Bove

simulations/missing-data-benchmarks/R/eval/helpers.R

@@ -1,3 +1,5 @@
+# get_trt_visit_num_ests ----
+
 # extract ATE estimates at each visit from mmrm fit
 get_mmrm_trt_visit_num_ests <- function(fit) {
   marginal_means <- emmeans::emmeans(
@@ -55,6 +57,8 @@ get_trt_visit_num_ests <- function(method, fit, dt, converged) {
   }
 }
 
+# get_convergence ----
+
 # extract convergence status from mmrm fit
 get_mmrm_convergence <- function(converged) {
   converged
@@ -88,6 +92,7 @@ get_convergence <- function(method, fit, converged) {
   }
 }
 
+# get_trt_visit_num_ses ----
 
 # extract standard errors for ATE estimates at each visit from mmrm fit
 get_mmrm_trt_visit_num_ses <- function(fit) {
@@ -146,6 +151,8 @@ get_trt_visit_num_ses <- function(method, fit, dt, converged) {
   }
 }
 
+# get_trt_visit_num_pvals ----
+
 # extract p-values for 2-sided hypothesis tests about ATEs at each visit from
 # mmrm fit
 get_mmrm_trt_visit_num_pvals <- function(fit) {
@@ -207,3 +214,112 @@ get_trt_visit_num_pvals <- function(method, fit, dt, converged) {
     rep(NA, 10) # hard code 10, the number of visits in all simulations
   }
 }
+
+# get_emmeans_output ----
+
+# compute the 95% CI from emmeans output
+get_95_ci <- function(emmeans_df) {
+  emmeans_df %>%
+    mutate(
+      lower = estimate - 1.96 * SE,
+      upper = estimate + 1.96 * SE,
+    )
+}
+
+format_emmeans_df <- function(emmeans_df) {
+  emmeans_df %>%
+    transmute(
+      visit_num = visit_num,
+      estimate = estimate,
+      stderr = SE,
+      df = df,
+      tvalue = t.ratio,
+      pvalue = p.value,
+      lower = lower,
+      upper = upper
+    )
+}
+
+# extract model fits output at each visit from mmrm fit
+get_mmrm_emmeans_output <- function(fit) {
+  # extract emmeans output
+  marginal_means <- emmeans(
+    fit,
+    spec = trt.vs.ctrl ~ trt | visit_num,
+    weights = "proportional"
+  )
+  emmeans_df <- as.data.frame(marginal_means$contrasts)
+
+  # compute lower and upper 95% CI
+  emmeans_df <- get_95_ci(emmeans_df)
+
+  # format to resemble SAS output
+  format_emmeans_df(emmeans_df)
+}
+
+# extract model fit outputs at each visit from glmmTMB fit
+get_glmm_emmeans_output <- function(fit) {
+  marginal_means <- emmeans(
+    fit,
+    spec = trt.vs.ctrl ~ trt | visit_num,
+    weights = "proportional"
+  )
+  emmeans_df <- as.data.frame(marginal_means$contrasts)
+
+  # compute lower and upper 95% CI
+  emmeans_df <- get_95_ci(emmeans_df)
+
+  # format to resemble SAS output
+  format_emmeans_df(emmeans_df)
+}
+
+# extract model fit outputs estimates at each visit from gls fit
+get_nlme_emmeans_output <- function(fit, data) {
+  marginal_means <- emmeans(
+    fit,
+    spec = trt.vs.ctrl ~ trt | visit_num,
+    weights = "proportional",
+    data = data,
+    mode = "df.error"
+  )
+  emmeans_df <- as.data.frame(marginal_means$contrasts)
+
+  # compute lower and upper 95% CI
+  emmeans_df <- get_95_ci(emmeans_df)
+
+  # format to resemble SAS output
+  format_emmeans_df(emmeans_df)
+}
+
+# extract model fit outputs estimates at each visit from PROC MIXED fit
+get_mixed_emmeans_like_output <- function(fit) {
+  fit %>% transmute(
+    visit_num = str_extract(contrast, "^([^:])+"),
+    estimate = Estimate,
+    stderr = StdErr,
+    df = DF,
+    tvalue = tValue,
+    pvalue = map2_dbl(tValue, DF, function(tvalue, df) {
+      2 * min(c(pt(tvalue, df), pt(tvalue, df, lower.tail = FALSE)))
+    }),
+    lower = Lower,
+    upper = Upper
+  )
+}
+
+# general function for emmeans like information
+get_emmeans_output <- function(method, fit, dt, converged) {
+  if (get_convergence(method, fit, converged)) {
+    if (str_detect(method, "mmrm")) {
+      get_mmrm_emmeans_output(fit)
+    } else if (str_detect(method, "glmmtmb")) {
+      get_glmm_emmeans_output(fit)
+    } else if (str_detect(method, "nlme")) {
+      get_nlme_emmeans_output(fit, dt)
+    } else if (str_detect(method, "proc_mixed")) {
+      get_mixed_emmeans_like_output(fit)
+    }
+  } else {
+    NA
+  }
+}

simulations/missing-data-benchmarks/tests/testthat.R

@@ -1,5 +1,6 @@
 # load required libraries
 library(simChef)
+library(checkmate)
 library(testthat)
 library(mmrm)
 library(glmmTMB)
@@ -11,6 +12,9 @@ library(tidyr)
 library(stringr)
 library(emmeans)
 
+# flag whether SAS tests should be run
+run_sas_tests <- FALSE
+
 # source the R scripts
 sim_functions_files <- list.files(
   c("R/dgp", "R/method", "R/eval", "R/viz"),

simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-bias.R

@@ -66,6 +66,9 @@ test_that("bias_fun computes DGP-specific bias", {
     bcva_change = eff_us$bcva_change,
     covar_type = "us"
   )
+
+  skip_if_not(run_sas_tests)
+
   proc_mixed_no_eff <- proc_mixed_wrapper_fun(
     participant = no_eff_us$participant,
     trt = no_eff_us$trt,

simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-convergence-rate.R

@@ -66,6 +66,9 @@ test_that("coconvergence_rate_rate_fun computes DGP-specific convergence rate",
     bcva_change = eff_us$bcva_change,
     covar_type = "us"
   )
+
+  skip_if_not(run_sas_tests)
+
   proc_mixed_no_eff <- proc_mixed_wrapper_fun(
     participant = no_eff_us$participant,
     trt = no_eff_us$trt,

simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-coverage.R

@@ -66,6 +66,9 @@ test_that("coverage_fun computes DGP-specific coverage", {
     bcva_change = eff_us$bcva_change,
     covar_type = "us"
   )
+
+  skip_if_not(run_sas_tests)
+
   proc_mixed_no_eff <- proc_mixed_wrapper_fun(
     participant = no_eff_us$participant,
     trt = no_eff_us$trt,

simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-helpers.R

@@ -1,3 +1,5 @@
+# get_trt_visit_num_ests ----
+
 test_that("get_mmrm_trt_visit_num_ests extracts interaction estimates", {
   # generate data
   set.seed(51235)
@@ -74,6 +76,7 @@ test_that("get_nlme_trt_visit_num_ests extracts interaction estimates", {
 })
 
 test_that("get_mixed_trt_visit_num_ests extracts interaction estimates", {
+  skip_if_not(run_sas_tests)
   # generate data
   set.seed(51235)
   dgp_output <- rct_dgp_fun(
@@ -98,6 +101,8 @@ test_that("get_mixed_trt_visit_num_ests extracts interaction estimates", {
   expect_equal(ests, c(0.5, 1.0, 1.5), tolerance = 0.1)
 })
 
+# get_convergence ----
+
 test_that("get_mmrm_convergence extracts convergence status", {
   # generate data
   set.seed(51235)
@@ -149,6 +154,7 @@ test_that("get_glmm_convergence extracts convergence status", {
 })
 
 test_that("get_mixed_convergence extracts convergence status", {
+  skip_if_not(run_sas_tests)
   # generate data
   set.seed(51235)
   dgp_output <- rct_dgp_fun(
@@ -173,6 +179,8 @@ test_that("get_mixed_convergence extracts convergence status", {
   expect_equal(conv_status, TRUE)
 })
 
+# get_trt_visit_num_ses ----
+
 test_that("get_mmrm_trt_visit_num_ses extracts standard errors", {
   # generate data
   set.seed(51235)
@@ -252,6 +260,7 @@ test_that("get_nlme_trt_visit_num_ses extracts standard errors", {
 })
 
 test_that("get_mixed_trt_visit_num_ses extracts standard errors", {
+  skip_if_not(run_sas_tests)
   # generate data
   set.seed(51235)
   dgp_output <- rct_dgp_fun(
@@ -277,6 +286,8 @@ test_that("get_mixed_trt_visit_num_ses extracts standard errors", {
   expect_equal(length(ses), 3)
 })
 
+# get_trt_visit_num_pvals ----
+
 test_that("get_mmrm_trt_visit_num_pvals extracts nominal p-values", {
   # generate data
   set.seed(51235)
@@ -356,6 +367,7 @@ test_that("get_nlme_trt_visit_num_pvals extracts nominal p-values", {
 })
 
 test_that("get_mixed_trt_visit_num_pvals extracts nominal p-values", {
+  skip_if_not(run_sas_tests)
   # generate data
   set.seed(51235)
   dgp_output <- rct_dgp_fun(
@@ -380,3 +392,118 @@ test_that("get_mixed_trt_visit_num_pvals extracts nominal p-values", {
   expect_equal(is.numeric(pvals), TRUE)
   expect_equal(length(pvals), 3)
 })
+
+# get_emmeans_output ----
+
+test_that("get_mmrm_emmeans_output extracts emmeans contrasts", {
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  mmrm_fit <- mmrm_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  # extract emmeans contrasts estimates
+  result <- get_mmrm_emmeans_output(mmrm_fit$fit)
+  expect_data_frame(result)
+  expect_identical(nrow(result), 3L)
+  expected_cols <- c("visit_num", "estimate", "stderr", "df", "tvalue", "pvalue", "lower", "upper")
+  expect_named(result, expected_cols)
+})
+
+test_that("get_glmm_emmeans_output extracts emmeans contrasts", {
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  glmmtmb_fit <- glmmtmb_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  # extract estimates
+  result <- get_glmm_emmeans_output(glmmtmb_fit$fit)
+  expect_data_frame(result)
+  expect_identical(nrow(result), 3L)
+  expected_cols <- c("visit_num", "estimate", "stderr", "df", "tvalue", "pvalue", "lower", "upper")
+  expect_named(result, expected_cols)
+})
+
+test_that("get_nlme_emmeans_output extracts emmeans contrasts", {
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  nlme_fit <- nlme_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  # extract estimates
+  result <- get_nlme_emmeans_output(nlme_fit$fit, nlme_fit$data)
+  expect_data_frame(result)
+  expect_identical(nrow(result), 3L)
+  expected_cols <- c("visit_num", "estimate", "stderr", "df", "tvalue", "pvalue", "lower", "upper")
+  expect_named(result, expected_cols)
+})
+
+test_that("get_mixed_emmeans_like_output extracts emmeans contrasts", {
+  skip_if_not(run_sas_tests)
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  proc_mixed_fit <- proc_mixed_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  # extract estimates
+  result <- get_mixed_emmeans_like_output(proc_mixed_fit$fit)
+  expect_data_frame(result)
+  expect_identical(nrow(result), 3L)
+  expected_cols <- c("visit_num", "estimate", "stderr", "df", "tvalue", "pvalue", "lower", "upper")
+  expect_named(result, expected_cols)
+})

simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-type-1-error.R

@@ -66,6 +66,7 @@ test_that("type_1_error_rate_fun computes DGP-specific type 1 error rates", {
     bcva_change = eff_us$bcva_change,
     covar_type = "us"
   )
+  skip_if_not(run_sas_tests)
   proc_mixed_no_eff <- proc_mixed_wrapper_fun(
     participant = no_eff_us$participant,
     trt = no_eff_us$trt,

simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-type-2-error.R

@@ -66,6 +66,9 @@ test_that("type_2_error_rate_fun computes DGP-specific type 2 error rates", {
     bcva_change = eff_us$bcva_change,
     covar_type = "us"
   )
+
+  skip_if_not(run_sas_tests)
+
   proc_mixed_no_eff <- proc_mixed_wrapper_fun(
     participant = no_eff_us$participant,
     trt = no_eff_us$trt,