add s_get_lsmeans to summarize the LS-means (#7)

.Rbuildignore

@@ -5,3 +5,4 @@
 ^LICENSE\.md$
 ^README\.Rmd$
 ^\.github$
+^data-raw$

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),

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)

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()
+}

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)
+}

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"
+  )
+}

inst/WORDLIST

@@ -1,4 +1,10 @@
 ADaM
+ANCOVA
 CDISC
+CMD
+MMRM
+PN
+TRT
 magrittr
 oversighting
+θ

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
+  )
+}