From 2567cced395f593b351c049c6abc7c04c40567c9 Mon Sep 17 00:00:00 2001 From: Kai Gu Date: Fri, 12 Jan 2024 10:07:08 +0800 Subject: [PATCH] add s_get_lsmeans to summarize the LS-means --- .Rbuildignore | 1 + DESCRIPTION | 7 +- NAMESPACE | 4 +- R/package.R | 31 ++++- R/pkg-methods.R | 35 ++++++ R/summarize-lsmeans.R | 152 ++++++++++++++++++++++++ inst/WORDLIST | 6 + man/s_get_lsmeans.Rd | 112 +++++++++++++++++ tests/spelling.R | 9 +- tests/testthat/test-pkg-methods.R | 20 ++++ tests/testthat/test-summarize-lsmeans.R | 111 +++++++++++++++++ 11 files changed, 480 insertions(+), 8 deletions(-) create mode 100644 R/pkg-methods.R create mode 100644 R/summarize-lsmeans.R create mode 100644 man/s_get_lsmeans.Rd create mode 100644 tests/testthat/test-pkg-methods.R create mode 100644 tests/testthat/test-summarize-lsmeans.R diff --git a/.Rbuildignore b/.Rbuildignore index 2fe6392..66be6e5 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,3 +5,4 @@ ^LICENSE\.md$ ^README\.Rmd$ ^\.github$ +^data-raw$ diff --git a/DESCRIPTION b/DESCRIPTION index 746273f..d4637df 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,12 +15,15 @@ Depends: R (>= 3.6) Imports: checkmate, + dplyr, + emmeans, lifecycle, magrittr, - methods, - stats + stats, + tibble Suggests: knitr, + mmrm, rmarkdown, spelling, testthat (>= 3.0.0), diff --git a/NAMESPACE b/NAMESPACE index 53e7570..f8342eb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,10 +1,12 @@ # Generated by roxygen2: do not edit by hand +S3method(print,s_lsmeans) export("%>%") export(rrPostProb) +export(s_get_lsmeans) import(checkmate) -import(methods) importFrom(lifecycle,deprecated) importFrom(magrittr,"%>%") +importFrom(stats,confint) importFrom(stats,pbeta) importFrom(stats,rbinom) diff --git a/R/package.R b/R/package.R index 81ecb10..8c96bd6 100644 --- a/R/package.R +++ b/R/package.R @@ -5,7 +5,34 @@ "_PACKAGE" #' @import checkmate -#' @import methods #' @importFrom lifecycle deprecated -#' @importFrom stats pbeta rbinom +#' @importFrom stats pbeta rbinom confint NULL + +.onLoad <- function(libname, pkgname) { + op <- options() + op.stabiot <- list( + stabiot.precision.default = tibble::tribble( + ~stat, ~extra, + "N", 0, + "MEAN", 1, + "SD", 2, + "MEDIAN", 1, + "MAX", 0, + "MIN", 0, + "Q1", 1, + "Q3", 1, + "EST", 1, + "SE", 2, + "DF", 0, + "CI", 1, + "RATIO", 1 + ) + ) + toset <- !(names(op.stabiot) %in% names(op)) + if (any(toset)) { + options(op.stabiot[toset]) + } + + invisible() +} diff --git a/R/pkg-methods.R b/R/pkg-methods.R new file mode 100644 index 0000000..9e085b2 --- /dev/null +++ b/R/pkg-methods.R @@ -0,0 +1,35 @@ +#' @describeIn s_get_lsmeans prints Least-squares Means summary. +#' @exportS3Method +#' @keywords internal +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 = "") + } + + cat("\n") + cat("Least-squares Means Estimates:\n") + print(x$lsm_est) + + 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) + } + + invisible(x) +} diff --git a/R/summarize-lsmeans.R b/R/summarize-lsmeans.R new file mode 100644 index 0000000..0cf3cbd --- /dev/null +++ b/R/summarize-lsmeans.R @@ -0,0 +1,152 @@ +# s_get_lsmeans ---- + +# s_get_lsmeans <- function(object, ...) { +# UseMethod("s_get_lsmeans") +# } + +#' Summarize Least-squares Means from Models +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' To estimate the Least-squares Means along with corresponding statistics and +#' confidence interval, and perform the hypothesis testing, using the `emmeans` +#' package. +#' +#' @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. +#' @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. +#' @param contrast (`logical`)\cr whether to perform the contrasts, pairwise +#' comparisons (required usage) and hypothesis testing or not, default is TRUE. +#' @param null (`numeric`)\cr null hypothesis value, it will also be referred to +#' the margin of superiority and non-inferiority design, default is 0. +#' @param conf.level (`numeric`)\cr significance level for the returned confidence +#' interval and hypothesis testing, default is 0.95. +#' @param alternative (`string`)\cr string specifying the alternative hypothesis, +#' must be one of "two.sided" (default), "greater" or "less". See the special +#' section below for more details. +#' @param ... other arguments to be passed to [emmeans::emmeans()]. +#' +#' @rdname s_get_lsmeans +#' @order 1 +#' +#' @return +#' An object of class `s_lsmeans` is a list contains several summary tables and +#' fit model as following components: +#' +#' - `model` information of fitted model. +#' - `lsm_est` estimate `tibble` table of Ls-means with corresponding statistics. +#' If `by` parameter is specified, all estimates for each `var` and `by` will +#' be presented. +#' - `lsm_contr` contrast `tibble` table with pairwise comparisons that contains +#' corresponding statistics, if the `contrast` is `TRUE`. +#' - `params` essential parameters. +#' +#' @note There is no any p value adjusted method involved. +#' +#' @details +#' For example, when the null hypothesis is θ equal to `0`, the `side` should be +#' `two.sided`, and `null` is defined as `0`. For the non-inferiority trial, +#' such as when the null hypothesis is θ inferiority to `-2`, the `side` should +#' be `greater`, and `null` is defined as `-2`. For the superiority trial, +#' when the null hypothesis is θ not superiority to `2`, the `side` should be +#' `greater`, and `null` is defined as `2`. +#' +#' @export +#' +#' @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::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)) %>% +#' lm(formula = FEV1 ~ FEV1_BL + RACE + SEX + ARMCD) +#' +#' s_get_lsmeans(fit2, "ARMCD") +s_get_lsmeans <- function(object, + var, + by = NULL, + contrast = TRUE, + null = 0, + conf.level = 0.95, + alternative = c("two.sided", "less", "greater"), + ...) { + assert_string(var) + assert_string(by, null.ok = TRUE) + assert_logical(contrast) + assert_number(null) + assert_number(conf.level, lower = 0, upper = 1) + alternative <- match.arg(alternative, c("two.sided", "less", "greater"), several.ok = FALSE) + + ems <- if (is.null(by)) { + emmeans::emmeans(object, var, ...) + } else { + emmeans::emmeans(object, var, by = by, ...) + } + + lsm_ci <- tibble::as_tibble( + data.frame(confint(ems, level = conf.level)) + ) + lsm_pval <- tibble::as_tibble( + data.frame(emmeans::test(ems)) + ) + lsm_est_res <- suppressMessages(dplyr::full_join(lsm_ci, lsm_pval)) %>% + dplyr::rename( + "estimate" = "emmean" + ) + + contr <- if (contrast) { + emmeans::contrast(ems, adjust = "none", method = "pairwise") + } + contr_ci <- tibble::as_tibble( + data.frame(confint(contr, level = conf.level)) + ) + side <- switch(alternative, + two.sided = 0, + less = -1, + greater = 1 + ) + contr_pval <- tibble::as_tibble( + data.frame(emmeans::test(contr, null = null, delta = 0, side = side, level = conf.level)) + ) + contr_est_res <- suppressMessages(dplyr::full_join(contr_ci, contr_pval)) + + structure( + list( + model = ems@model.info, + lsm_est = lsm_est_res, + lsm_contr = contr_est_res, + params = list( + var = var, + by = by, + contrast = contrast, + null = null, + conf.level = conf.level, + alternative = alternative + ) + ), + class = "s_lsmeans" + ) +} diff --git a/inst/WORDLIST b/inst/WORDLIST index fbb138e..f2971c8 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,4 +1,10 @@ ADaM +ANCOVA CDISC +CMD +MMRM +PN +TRT magrittr oversighting +θ diff --git a/man/s_get_lsmeans.Rd b/man/s_get_lsmeans.Rd new file mode 100644 index 0000000..971b095 --- /dev/null +++ b/man/s_get_lsmeans.Rd @@ -0,0 +1,112 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/summarize-lsmeans.R, R/pkg-methods.R +\name{s_get_lsmeans} +\alias{s_get_lsmeans} +\alias{print.s_lsmeans} +\title{Summarize Least-squares Means from Models} +\usage{ +s_get_lsmeans( + object, + var, + by = NULL, + contrast = TRUE, + null = 0, + conf.level = 0.95, + alternative = c("two.sided", "less", "greater"), + ... +) + +\method{print}{s_lsmeans}(x, ...) +} +\arguments{ +\item{object}{(\verb{fitted model})\cr any fitted model that can be accepted by +\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.} + +\item{by}{(\code{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.} + +\item{contrast}{(\code{logical})\cr whether to perform the contrasts, pairwise +comparisons (required usage) and hypothesis testing or not, default is TRUE.} + +\item{null}{(\code{numeric})\cr null hypothesis value, it will also be referred to +the margin of superiority and non-inferiority design, default is 0.} + +\item{conf.level}{(\code{numeric})\cr significance level for the returned confidence +interval and hypothesis testing, default is 0.95.} + +\item{alternative}{(\code{string})\cr string specifying the alternative hypothesis, +must be one of "two.sided" (default), "greater" or "less". See the special +section below for more details.} + +\item{...}{other arguments to be passed to \code{\link[emmeans:emmeans]{emmeans::emmeans()}}.} +} +\value{ +An object of class \code{s_lsmeans} is a list contains several summary tables and +fit model as following components: +\itemize{ +\item \code{model} information of fitted model. +\item \code{lsm_est} estimate \code{tibble} table of Ls-means with corresponding statistics. +If \code{by} parameter is specified, all estimates for each \code{var} and \code{by} will +be presented. +\item \code{lsm_contr} contrast \code{tibble} table with pairwise comparisons that contains +corresponding statistics, if the \code{contrast} is \code{TRUE}. +\item \code{params} essential parameters. +} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +To estimate the Least-squares Means along with corresponding statistics and +confidence interval, and perform the hypothesis testing, using the \code{emmeans} +package. +} +\details{ +For example, when the null hypothesis is θ equal to \code{0}, the \code{side} should be +\code{two.sided}, and \code{null} is defined as \code{0}. For the non-inferiority trial, +such as when the null hypothesis is θ inferiority to \code{-2}, the \code{side} should +be \code{greater}, and \code{null} is defined as \code{-2}. For the superiority trial, +when the null hypothesis is θ not superiority to \code{2}, the \code{side} should be +\code{greater}, and \code{null} is defined as \code{2}. +} +\section{Functions}{ +\itemize{ +\item \code{print(s_lsmeans)}: prints Least-squares Means summary. + +}} +\note{ +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::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)) \%>\% + lm(formula = FEV1 ~ FEV1_BL + RACE + SEX + ARMCD) + +s_get_lsmeans(fit2, "ARMCD") +} +\keyword{internal} 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 new file mode 100644 index 0000000..ff69bbe --- /dev/null +++ b/tests/testthat/test-pkg-methods.R @@ -0,0 +1,20 @@ +test_that("print.s_lsmeans works as expected for mmrm", { + data(fev_data, package = "mmrm") + fit <- mmrm::mmrm( + formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID), + reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear", + data = fev_data + ) + 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) + + 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) + + 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) +}) diff --git a/tests/testthat/test-summarize-lsmeans.R b/tests/testthat/test-summarize-lsmeans.R new file mode 100644 index 0000000..a4107d3 --- /dev/null +++ b/tests/testthat/test-summarize-lsmeans.R @@ -0,0 +1,111 @@ +test_that("s_get_lsmeans works as expected by visit with default arguments", { + data(fev_data, package = "mmrm") + fit <- mmrm::mmrm( + formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID), + reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear", + data = fev_data + ) + res <- s_get_lsmeans(fit, "ARMCD", by = "AVISIT") + + expect_class(res, "s_lsmeans") + expect_identical(dim(res$lsm_est), c(8L, 9L)) + expect_identical(dim(res$lsm_contr), c(4L, 9L)) + expect_equal( + res$lsm_est[1:2, ], + tibble::tibble( + ARMCD = factor(c("PBO", "TRT")), + AVISIT = factor(c("VIS1", "VIS1"), levels = c("VIS1", "VIS2", "VIS3", "VIS4")), + estimate = c(33.33186, 37.10609), + SE = c(0.7612042, 0.7674844), + df = c(148.1457, 143.1765), + lower.CL = c(31.82764, 35.58903), + upper.CL = c(34.83608, 38.62316), + t.ratio = c(43.78833, 48.34768), + p.value = c(1.165649e-86, 1.444629e-90) + ), + tolerance = 0.0001 + ) + expect_equal( + res$lsm_contr[1:2, ], + tibble::tibble( + contrast = factor(c("PBO - TRT", "PBO - TRT")), + AVISIT = factor(c("VIS1", "VIS2"), levels = c("VIS1", "VIS2", "VIS3", "VIS4")), + 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), + p.value = c(6.417790e-04, 2.840654e-05) + ), + tolerance = 0.0001 + ) +}) + +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", + data = fev_data + ) + res <- s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = 2, alternative = "greater") + + expect_identical(dim(res$lsm_contr), c(4L, 10L)) + expect_equal( + res$lsm_contr[1:2, ], + tibble::tibble( + contrast = factor(c("TRT - PBO", "TRT - PBO")), + AVISIT = factor(c("VIS1", "VIS2"), levels = c("VIS1", "VIS2", "VIS3", "VIS4")), + estimate = c(3.774230, 3.732304), + SE = c(1.0817635, 0.8633334), + df = c(145.5520, 145.2771), + lower.CL = c(1.636236, 2.025988), + upper.CL = c(5.912224, 5.438620), + null = c(2, 2), + t.ratio = c(1.640127, 2.006529), + p.value = c(0.05156904, 0.02332779) + ), + tolerance = 0.0001 + ) +}) + +test_that("s_get_lsmeans works as expected for ancova", { + data(fev_data, package = "mmrm") + fit <- fev_data %>% + dplyr::filter(VISITN == 4 & !is.na(FEV1)) %>% + lm(formula = FEV1 ~ FEV1_BL + RACE + SEX + ARMCD) + res <- s_get_lsmeans(fit, "ARMCD") + + expect_identical(dim(res$lsm_est), c(2L, 8L)) + expect_identical(dim(res$lsm_contr), c(1L, 8L)) + expect_equal( + res$lsm_est, + tibble::tibble( + ARMCD = factor(c("PBO", "TRT")), + estimate = c(48.63345, 52.53288), + SE = c(1.204272, 1.180020), + df = c(128, 128), + lower.CL = c(46.25060, 50.19801), + upper.CL = c(51.01631, 54.86775), + t.ratio = c(40.38412, 44.51864), + p.value = c(1.071522e-74, 9.320133e-80) + ), + tolerance = 0.0001 + ) + expect_equal( + res$lsm_contr, + tibble::tibble( + contrast = c("PBO - TRT"), + 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), + p.value = c(0.02235472) + ), + tolerance = 0.0001 + ) +})