From 44017f7beb7adb7c1a3f9f0ecb05b7e2edc77ac2 Mon Sep 17 00:00:00 2001 From: Kai Gu Date: Wed, 24 Apr 2024 15:59:20 +0800 Subject: [PATCH] update LSmean print method and tests (#23) --- NAMESPACE | 1 + R/package.R | 4 +- R/pkg-methods.R | 93 +++++++++++++++++++------ R/summarize-lsmeans.R | 13 ++-- man/s_get_lsmeans.Rd | 10 +-- tests/spelling.R | 9 ++- tests/testthat/test-pkg-methods.R | 10 ++- tests/testthat/test-summarize-lsmeans.R | 21 +++--- 8 files changed, 103 insertions(+), 58 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 16ddbdd..f7ffbaa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,6 +19,7 @@ export(s_propci) import(checkmate) importFrom(dplyr,add_count) importFrom(dplyr,arrange) +importFrom(dplyr,bind_rows) importFrom(dplyr,case_when) importFrom(dplyr,count) importFrom(dplyr,distinct) diff --git a/R/package.R b/R/package.R index b37d535..32ce4bd 100644 --- a/R/package.R +++ b/R/package.R @@ -8,8 +8,8 @@ #' @importFrom magrittr set_colnames set_rownames #' @importFrom lifecycle deprecated #' @importFrom stats pbeta rbinom confint as.formula setNames coef quantile -#' @importFrom dplyr add_count arrange case_when count distinct filter full_join -#' group_by left_join mutate row_number rowwise select summarise ungroup +#' @importFrom dplyr add_count arrange bind_rows case_when count distinct filter +#' full_join group_by left_join mutate row_number rowwise select summarise ungroup #' @importFrom rlang sym := .data #' @importFrom rtables analyze basic_table build_table in_rows rcell non_ref_rcell #' split_cols_by diff --git a/R/pkg-methods.R b/R/pkg-methods.R index 8f65410..b6a822c 100644 --- a/R/pkg-methods.R +++ b/R/pkg-methods.R @@ -5,32 +5,83 @@ print.s_lsmeans <- function(x, ...) { cat("Model Call: ", append = FALSE) print(x$model$call) cat("\n") - varc <- paste0(x$model$xlev[[x$params$var]], collapse = ", ") - cat("Predictor/Treatment: ", x$params$var, " (", varc, ")\n", sep = "") - if (!is.null(x$params$by)) { - byc <- paste0(x$model$xlev[[x$params$by]], collapse = ", ") - cat("Group by: ", x$params$by, " (", byc, ")\n", sep = "") + + grp_var <- x$params$var + by_var <- x$params$by + df <- x$data + ref_col <- levels(x$lsm_est[[grp_var]])[1] + side <- if (x$params$alternative == "greater" & x$params$null >= 0) { + "p-value (Superiority Test)" + } else if (x$params$alternative == "greater" & x$params$null < 0) { + "p-value (Non-inferiority Test)" + } else if (x$params$alternative == "less" & x$params$null >= 0) { + "p-value (Non-Superiority Test)" + } else if (x$params$alternative == "less" & x$params$null < 0) { + "p-value (inferiority Test)" + } else { + "p-value" } - cat("\n") - cat("Least-squares Means Estimates:\n") - print(x$lsm_est) + df <- bind_rows( + x$lsm_est, + x$lsm_contr + ) %>% + rowwise() %>% + mutate( + !!sym(grp_var) := case_when( + is.na(!!sym(grp_var)) ~ strsplit(as.character(.data$contrast), " - ")[[1]][1], + TRUE ~ !!sym(grp_var) + ) + ) - if (x$params$contrast) { - cat("\n") - cat("Contrast Estimates of Least-squares Means:\n") - if (x$params$alternative == "two.sided") { - cat(sprintf("Null hypothesis is \u03b8 equal to %s.\n", x$params$null)) - } else if (x$params$alternative == "greater" & x$params$null >= 0) { - cat(sprintf("Null hypothesis is \u03b8 non-superiority to %s.\n", x$params$null)) - } else if (x$params$alternative == "greater" & x$params$null < 0) { - cat(sprintf("Null hypothesis is \u03b8 inferiority to %s.\n", x$params$null)) - } else if (x$params$alternative == "less" & x$params$null >= 0) { - cat(sprintf("Null hypothesis is \u03b8 superiority to %s.\n", x$params$null)) - } - print(x$lsm_contr) + a_lsmeam_func <- function(df, .var, .in_ref_col, est, contr, vst) { + curgrp <- df[[.var]][1] + df_est <- df %>% + filter(!!sym(grp_var) == curgrp & is.na(.data$contrast)) + df_contr <- df %>% + filter(!!sym(grp_var) == curgrp & .data$contrast == paste(curgrp, "-", ref_col)) + in_rows( + rcell(unlist(df_est[, c("estimate", "SE"), drop = TRUE]), format = "xx.xx (xx.xx)"), + rcell(unlist(df_est[, c("lower.CL", "upper.CL"), drop = TRUE]), + format = "(xx.xx, xx.xx)", + indent_mod = 1L + ), + non_ref_rcell( + unlist(df_contr[, c("estimate", "SE"), drop = TRUE]), + .in_ref_col, + format = "xx.xx (xx.xx)" + ), + non_ref_rcell( + unlist(df_contr[, c("lower.CL", "upper.CL"), drop = TRUE]), + .in_ref_col, + format = "(xx.xx, xx.xx)", + indent_mod = 1L + ), + non_ref_rcell( + df_contr[, "p.value", drop = TRUE], + .in_ref_col, + format = "x.xxxx | (<0.0001)", + indent_mod = 1L + ), + .names = c( + "LS Mean (SE)", "95% CI", + "Difference LS Mean (SE)", "95% CI", + side + ) + ) } + tbl <- basic_table() %>% + split_cols_by(grp_var, ref_group = ref_col) + if (!is.null(by_var)) { + tbl <- tbl %>% rtables::split_rows_by(by_var, label_pos = "topleft") + } + result <- tbl %>% + analyze(grp_var, a_lsmeam_func) %>% + rtables::append_topleft(" Statistics") %>% + build_table(df = df) + print(result) + invisible(x) } diff --git a/R/summarize-lsmeans.R b/R/summarize-lsmeans.R index 45b68ea..3eed71c 100644 --- a/R/summarize-lsmeans.R +++ b/R/summarize-lsmeans.R @@ -15,7 +15,8 @@ #' @param object (`fitted model`)\cr any fitted model that can be accepted by #' `emmeans::emmeans()`, such as object from `lm` for ANCOVA and `mmrm` for MMRM. #' @param var (`string`)\cr string specifying the name of the predictor, such as -#' it would be treatment group variable (TRT01PN) in the clinical trials. +#' it would be treatment group variable (TRT01PN) in the clinical trials. Default +#' the first level is defined as reference group. #' @param by (`string`)\cr string specifying the name of the grouped variable. #' The estimates and contrasts will be evaluated separately for each elements #' from the grouped variable. @@ -74,25 +75,20 @@ #' @examples #' # refactor the level order: #' data(fev_data, package = "mmrm") -#' fev_data$ARMCD <- factor(fev_data$ARMCD, level = c("TRT", "PBO")) #' -#' # fit mmrm model: +#' # fit MMRM model: #' fit <- mmrm::mmrm( #' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID), #' reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear", #' data = fev_data #' ) #' -#' # estimate ARMCD within whole visits: -#' s_get_lsmeans(fit, "ARMCD") -#' #' # estimate ARMCD by visits: #' s_get_lsmeans(fit, "ARMCD", by = "AVISIT") #' #' # estimate ARMCD by visits for superiority testing with null hypothesis of 2: #' s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = 2, alternative = "greater") #' -#' #' # fit ANCOVA model: #' fit2 <- fev_data |> #' dplyr::filter(VISITN == 4 & !is.na(FEV1)) |> @@ -132,7 +128,7 @@ s_get_lsmeans <- function(object, ) contr <- if (contrast) { - emmeans::contrast(ems, adjust = "none", method = "pairwise") + emmeans::contrast(ems, adjust = "none", method = "revpairwise") } contr_ci <- tibble::as_tibble( data.frame(confint(contr, level = conf.level)) @@ -150,6 +146,7 @@ s_get_lsmeans <- function(object, structure( list( model = ems@model.info, + data = object$data, lsm_est = lsm_est_res, lsm_contr = contr_est_res, params = list( diff --git a/man/s_get_lsmeans.Rd b/man/s_get_lsmeans.Rd index 5e27294..14c7e66 100644 --- a/man/s_get_lsmeans.Rd +++ b/man/s_get_lsmeans.Rd @@ -23,7 +23,8 @@ s_get_lsmeans( \code{emmeans::emmeans()}, such as object from \code{lm} for ANCOVA and \code{mmrm} for MMRM.} \item{var}{(\code{string})\cr string specifying the name of the predictor, such as -it would be treatment group variable (TRT01PN) in the clinical trials.} +it would be treatment group variable (TRT01PN) in the clinical trials. Default +the first level is defined as reference group.} \item{by}{(\code{string})\cr string specifying the name of the grouped variable. The estimates and contrasts will be evaluated separately for each elements @@ -83,25 +84,20 @@ There is no any p value adjusted method involved. \examples{ # refactor the level order: data(fev_data, package = "mmrm") -fev_data$ARMCD <- factor(fev_data$ARMCD, level = c("TRT", "PBO")) -# fit mmrm model: +# fit MMRM model: fit <- mmrm::mmrm( formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID), reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear", data = fev_data ) -# estimate ARMCD within whole visits: -s_get_lsmeans(fit, "ARMCD") - # estimate ARMCD by visits: s_get_lsmeans(fit, "ARMCD", by = "AVISIT") # estimate ARMCD by visits for superiority testing with null hypothesis of 2: s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = 2, alternative = "greater") - # fit ANCOVA model: fit2 <- fev_data |> dplyr::filter(VISITN == 4 & !is.na(FEV1)) |> diff --git a/tests/spelling.R b/tests/spelling.R index 6713838..13f77d9 100644 --- a/tests/spelling.R +++ b/tests/spelling.R @@ -1,3 +1,6 @@ -if(requireNamespace('spelling', quietly = TRUE)) - spelling::spell_check_test(vignettes = TRUE, error = FALSE, - skip_on_cran = TRUE) +if (requireNamespace("spelling", quietly = TRUE)) { + spelling::spell_check_test( + vignettes = TRUE, error = FALSE, + skip_on_cran = TRUE + ) +} diff --git a/tests/testthat/test-pkg-methods.R b/tests/testthat/test-pkg-methods.R index d93c4df..3a6bfed 100644 --- a/tests/testthat/test-pkg-methods.R +++ b/tests/testthat/test-pkg-methods.R @@ -7,16 +7,14 @@ test_that("print.s_lsmeans works as expected for mmrm", { ) res <- capture_output(print(s_get_lsmeans(fit, "ARMCD", by = "AVISIT"))) expect_match(res, "Model Call: mmrm::mmrm(formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT", fixed = TRUE) - expect_match(res, "Predictor/Treatment: ARMCD (PBO, TRT)", fixed = TRUE) - expect_match(res, "Group by: AVISIT (VIS1, VIS2, VIS3, VIS4)", fixed = TRUE) - expect_match(res, "Least-squares Means Estimates", fixed = TRUE) - expect_match(res, "Null hypothesis is θ equal to 0", fixed = TRUE) + expect_match(res, "LS Mean (SE)", fixed = TRUE) + expect_match(res, "Difference LS Mean (SE)", fixed = TRUE) res2 <- capture_output(print(s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = 2, alternative = "greater"))) - expect_match(res2, "Null hypothesis is θ non-superiority to 2", fixed = TRUE) + expect_match(res2, "p-value (Superiority Test)", fixed = TRUE) res3 <- capture_output(print(s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = -2, alternative = "greater"))) - expect_match(res3, "Null hypothesis is θ inferiority to -2", fixed = TRUE) + expect_match(res3, "p-value (Non-inferiority Test)", fixed = TRUE) }) test_that("print.prop_ci works as expected", { diff --git a/tests/testthat/test-summarize-lsmeans.R b/tests/testthat/test-summarize-lsmeans.R index a4107d3..92ed704 100644 --- a/tests/testthat/test-summarize-lsmeans.R +++ b/tests/testthat/test-summarize-lsmeans.R @@ -28,14 +28,14 @@ test_that("s_get_lsmeans works as expected by visit with default arguments", { expect_equal( res$lsm_contr[1:2, ], tibble::tibble( - contrast = factor(c("PBO - TRT", "PBO - TRT")), + contrast = factor(c("TRT - PBO", "TRT - PBO")), AVISIT = factor(c("VIS1", "VIS2"), levels = c("VIS1", "VIS2", "VIS3", "VIS4")), - estimate = c(-3.774230, -3.732304), + estimate = c(3.774230, 3.732304), SE = c(1.0817635, 0.8633334), df = c(145.5520, 145.2771), - lower.CL = c(-5.912224, -5.438620), - upper.CL = c(-1.636236, -2.025988), - t.ratio = c(-3.488960, -4.323131), + lower.CL = c(1.636236, 2.025988), + upper.CL = c(5.912224, 5.438620), + t.ratio = c(3.488960, 4.323131), p.value = c(6.417790e-04, 2.840654e-05) ), tolerance = 0.0001 @@ -44,7 +44,6 @@ test_that("s_get_lsmeans works as expected by visit with default arguments", { test_that("s_get_lsmeans works as expected for superiority testing with null hypothesis of 2", { data(fev_data, package = "mmrm") - fev_data$ARMCD <- factor(fev_data$ARMCD, level = c("TRT", "PBO")) fit <- mmrm::mmrm( formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID), reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear", @@ -97,13 +96,13 @@ test_that("s_get_lsmeans works as expected for ancova", { expect_equal( res$lsm_contr, tibble::tibble( - contrast = c("PBO - TRT"), - estimate = c(-3.89943), + contrast = c("TRT - PBO"), + estimate = c(3.89943), SE = c(1.686347), df = c(128), - lower.CL = c(-7.236155), - upper.CL = c(-0.562705), - t.ratio = c(-2.312354), + lower.CL = c(0.562705), + upper.CL = c(7.236155), + t.ratio = c(2.312354), p.value = c(0.02235472) ), tolerance = 0.0001