update anova type II and III (#391)

DESCRIPTION

@@ -57,6 +57,7 @@ Imports:
     TMB (>= 1.9.1),
     utils
 Suggests:
+    car (>= 3.1.2),
     cli,
     clubSandwich,
     clusterGeneration,
@@ -119,4 +120,5 @@ Collate:
     'tmb.R'
     'interop-emmeans.R'
     'interop-parsnip.R'
+    'interop-car.R'
     'zzz.R'

NEWS.md

@@ -1,9 +1,14 @@
 # mmrm 0.3.7.9004
 
+### New Features
+
+- `Anova` is implemented for `mmrm` models and available upon loading the `car` package. It supports type II and III hypothesis testing.
+
 ### Miscellaneous
 
 - In documentation of `mmrm_control()`, the allowed vcov definition is corrected to "Empirical-Jackknife" (CR3), and "Empirical-Bias-Reduced" (CR2).
 - Fix a compiler warning related to missing format specification in error message function call.
+- If an empty contrast matrix is provided to `df_md`, it will return statistics with `NA` values.
 
 # mmrm 0.3.6
 

R/interop-car.R

@@ -0,0 +1,114 @@
+#' Register `mmrm` For Use With `car::Anova`
+#'
+#' @inheritParams base::requireNamespace
+#' @return A logical value indicating whether registration was successful.
+#'
+#' @keywords internal
+car_add_mmrm <- function(quietly = FALSE) {
+  if (!requireNamespace("car", quietly = quietly)) {
+    return(FALSE)
+  }
+  envir <- asNamespace("mmrm")
+  h_register_s3("car", "Anova", "mmrm", envir)
+  TRUE
+}
+
+
+#' Obtain Contrast for Specified Effect
+#'
+#' This is support function to obtain contrast matrix for type II/III testing.
+#'
+#' @param object (`mmrm`)\cr the fitted MMRM.
+#' @param effect (`string`) the name of the effect.
+#' @param type (`string`) type of test, "II", "III", '2', or '3'.
+#' @param tol (`numeric`) threshold blow which values are treated as 0.
+#'
+#' @return A `matrix` of the contrast.
+#'
+#' @keywords internal
+h_get_contrast <- function(object, effect, type = c("II", "III", "2", "3"), tol = sqrt(.Machine$double.eps)) {
+  assert_class(object, "mmrm")
+  assert_string(effect)
+  assert_double(tol, finite = TRUE, len = 1L)
+  type <- match.arg(type)
+  mx <- component(object, "x_matrix")
+  asg <- attr(mx, "assign")
+  formula <- object$formula_parts$model_formula
+  tms <- terms(formula)
+  fcts <- attr(tms, "factors")[-1L, , drop = FALSE] # Discard the response.
+  ods <- attr(tms, "order")
+  assert_subset(effect, colnames(fcts))
+  idx <- which(effect == colnames(fcts))
+  cols <- which(asg == idx)
+  data <- component(object, "full_frame")
+  var_numeric <- vapply(data, is.numeric, FUN.VALUE = TRUE)
+  coef_rows <- length(cols)
+  l_mx <- matrix(0, nrow = coef_rows, ncol = length(asg))
+  if (coef_rows == 0L) {
+    return(l_mx)
+  }
+  l_mx[, cols] <- diag(rep(1, coef_rows))
+  for (i in setdiff(seq_len(ncol(fcts)), idx)) {
+    x1 <- mx[, cols, drop = FALSE]
+    additional_vars <- names(which(fcts[, i] > fcts[, idx]))
+    additional_numeric <- any(var_numeric[additional_vars])
+    current_col <- which(asg == i)
+    if (ods[i] >= ods[idx] && all(fcts[, i] >= fcts[, idx]) && !additional_numeric) {
+      sub_mat <- switch(type,
+        "2" = ,
+        "II" = {
+          x0 <- mx[, -c(cols, current_col), drop = FALSE]
+          x2 <- mx[, current_col, drop = FALSE]
+          m <- diag(rep(1, nrow(x0))) - x0 %*% solve(t(x0) %*% x0) %*% t(x0)
+          solve(t(x1) %*% m %*% x1) %*% t(x1) %*% m %*% x2
+        },
+        "3" = ,
+        "III" = {
+          additional_levels <- vapply(data[additional_vars], function(x) nlevels(x), FUN.VALUE = 1L)
+          t_levels <- prod(additional_levels)
+          l_mx[, cols] / t_levels
+        }
+      )
+      l_mx[, current_col] <- sub_mat
+    }
+  }
+  l_mx[abs(l_mx) < tol] <- 0
+  l_mx
+}
+
+#' Conduct type II/III hypothesis testing on the MMRM fit results.
+#'
+#' @param mod (`mmrm`)\cr the fitted MMRM.
+#' @param ... not used.
+#' @inheritParams h_get_contrast
+#'
+#' @details
+#' `Anova` will return `anova` object with one row per variable and columns
+#' `Num Df`(numerator degrees of freedom), `Denom Df`(denominator degrees of freedom),
+#' `F Statistic` and `Pr(>=F)`.
+#'
+#' @keywords internal
+# Please do not load `car` and then create the documentation. The Rd file will be different.
+Anova.mmrm <- function(mod, type = c("II", "III", "2", "3"), tol = sqrt(.Machine$double.eps), ...) { # nolint
+  assert_double(tol, finite = TRUE, len = 1L)
+  type <- match.arg(type)
+  vars <- colnames(attr(terms(mod$formula_parts$model_formula), "factors"))
+  ret <- lapply(
+    vars,
+    function(x) df_md(mod, h_get_contrast(mod, x, type, tol))
+  )
+  ret_df <- do.call(rbind.data.frame, ret)
+  row.names(ret_df) <- vars
+  colnames(ret_df) <- c("Num Df", "Denom Df", "F Statistic", "Pr(>=F)")
+  class(ret_df) <- c("anova", "data.frame")
+  attr(ret_df, "heading") <- sprintf(
+    "Analysis of Fixed Effect Table (Type %s F tests)",
+    switch(type,
+      "2" = ,
+      "II" = "II",
+      "3" = ,
+      "III" = "III"
+    )
+  )
+  ret_df
+}

R/testing.R

@@ -64,6 +64,16 @@ df_md <- function(object, contrast) {
     contrast <- matrix(contrast, ncol = length(contrast))
   }
   assert_matrix(contrast, ncols = length(component(object, "beta_est")))
+  if (nrow(contrast) == 0) {
+    return(
+      list(
+        num_df = 0,
+        denom_df = NA_real_,
+        f_stat = NA_real_,
+        p_val = NA_real_
+      )
+    )
+  }
   switch(object$method,
     "Satterthwaite" = h_df_md_sat(object, contrast),
     "Kenward-Roger" = h_df_md_kr(object, contrast),

R/utils.R

@@ -347,3 +347,28 @@ h_valid_formula <- function(formula) {
     stop("`.` is not allowed in mmrm models!")
   }
 }
+
+#' Register S3 Method
+#' Register S3 method to a generic.
+#'
+#' @param pkg (`string`) name of the package name.
+#' @param generic (`string`) name of the generic.
+#' @param class (`string`) class name the function want to dispatch.
+#' @param envir (`environment`) the location the method is defined.
+#'
+#' @details This function is adapted from `emmeans:::register_s3_method()`.
+#'
+#' @keywords internal
+h_register_s3 <- function(pkg, generic, class, envir = parent.frame()) {
+  assert_string(pkg)
+  assert_string(generic)
+  assert_string(class)
+  assert_environment(envir)
+  fun <- get(paste0(generic, ".", class), envir = envir)
+  if (isNamespaceLoaded(pkg)) {
+    registerS3method(generic, class, fun, envir = asNamespace(pkg))
+  }
+  setHook(packageEvent(pkg, "onLoad"), function(...) {
+    registerS3method(generic, class, fun, envir = asNamespace(pkg))
+  })
+}

R/zzz.R

@@ -15,6 +15,11 @@
     callback = parsnip_add_mmrm,
     message = emit_tidymodels_register_msg
   )
+  register_on_load(
+    "car", c("3.1.2", NA),
+    callback = car_add_mmrm,
+    message = "mmrm() registered as car::Anova extension"
+  )
 }
 
 #' Helper Function for Registering Functionality With Suggests Packages

_pkgdown.yml

@@ -33,9 +33,10 @@ articles:
     navbar: Covariance Structures
     contents:
       - covariance
-  - title: Degrees of Freedom
-    navbar: Degrees of Freedom
+  - title: Degrees of Freedom and Testing
+    navbar: Degrees of Freedom and Testing
     contents:
+      - hypothesis_testing
       - between_within
       - kenward
       - satterthwaite

design/anova/anova.Rmd

@@ -533,3 +533,38 @@ df_md(fit, a)
 
 get_l_matrix(fit, "FEV1_BL", "2")
 ```
+
+## SAS Results for Integration Test
+
+```{r}
+df2sd(fev_data, "fev")
+sas_code <- "
+ods output Tests1 = Tests1 Tests2 = Tests2 Tests3 = Tests3;
+PROC MIXED DATA = fev cl method=reml;
+  CLASS AVISIT(ref = last) ARMCD(ref = 'PBO') USUBJID;
+  MODEL FEV1 = ARMCD FEV1_BL ARMCD * FEV1_BL / ddfm=satterthwaite htype=1,2,3 solution chisq e3;
+  REPEATED AVISIT / subject=USUBJID type=ar(1);
+RUN;
+"
+res <- run_sas(sas_code)
+df2 <- sd2df("tests2")
+df3 <- sd2df("tests3")
+write.csv(df2, "design/anova/test2_1.csv")
+write.csv(df3, "design/anova/test3_1.csv")
+```
+
+```{r}
+sas_code <- "
+ods output Tests1 = Tests1 Tests2 = Tests2 Tests3 = Tests3;
+PROC MIXED DATA = fev cl method=reml;
+  CLASS AVISIT(ref = last) SEX(ref = 'Male') ARMCD(ref = 'PBO') USUBJID;
+  MODEL FEV1 = ARMCD SEX ARMCD * SEX ARMCD * FEV1_BL / ddfm=satterthwaite htype=1,2,3 solution chisq e3;
+  REPEATED AVISIT / subject=USUBJID type=ar(1);
+RUN;
+"
+res <- run_sas(sas_code)
+df2 <- sd2df("tests2")
+df3 <- sd2df("tests3")
+write.csv(df2, "design/anova/test2_2.csv")
+write.csv(df3, "design/anova/test3_2.csv")
+```

design/anova/test2_1.csv

@@ -0,0 +1,4 @@
+"","Effect","NumDF","DenDF","ChiSq","FValue","ProbChiSq","ProbF"
+"1","ARMCD",1,199.582314485653,0.03099660681111,0.03099660681111,0.86024795837881,0.86042646803505
+"2","FEV1_BL",1,196.41571750149,11.5301174147263,11.5301174147263,0.00068477688975,0.00082860861048
+"3","FEV1_BL*ARMCD",1,198.163532037944,0.7644368094469,0.7644368094469,0.38194359117097,0.38300197080559

design/anova/test2_2.csv

@@ -0,0 +1,5 @@
+"","Effect","NumDF","DenDF","ChiSq","FValue","ProbChiSq","ProbF"
+"1","ARMCD",1,196.700026021882,0.03093826913613,0.03093826913613,0.86037817970071,0.86055912789767
+"2","SEX",1,186.404412481085,0.10489707559679,0.10489707559679,0.74603026404544,0.74639336822641
+"3","SEX*ARMCD",1,186.11901853766,0.46643535342468,0.46643535342468,0.49463165974308,0.49548089306181
+"4","FEV1_BL*ARMCD",2,194.942672335072,12.2296433517045,6.11482167585228,0.00220986980743,0.0026566780461

design/anova/test3_1.csv

@@ -0,0 +1,4 @@
+"","Effect","NumDF","DenDF","ChiSq","FValue","ProbChiSq","ProbF"
+"1","ARMCD",1,199.582314485653,0.03099660681111,0.03099660681111,0.86024795837881,0.86042646803505
+"2","FEV1_BL",1,198.163532037936,11.0435094615508,11.0435094615508,0.00088998227841,0.00105994994938
+"3","FEV1_BL*ARMCD",1,198.163532037944,0.7644368094469,0.7644368094469,0.38194359117097,0.38300197080559

design/anova/test3_2.csv

@@ -0,0 +1,5 @@
+"","Effect","NumDF","DenDF","ChiSq","FValue","ProbChiSq","ProbF"
+"1","ARMCD",1,197.423866403271,0.04875135180324,0.04875135180324,0.82525043506483,0.82547861285648
+"2","SEX",1,186.119018537661,0.09277653315012,0.09277653315012,0.76067660529697,0.76101691812013
+"3","SEX*ARMCD",1,186.11901853766,0.46643535342468,0.46643535342468,0.49463165974308,0.49548089306181
+"4","FEV1_BL*ARMCD",2,194.942672335072,12.2296433517045,6.11482167585228,0.00220986980743,0.0026566780461

inst/REFERENCES.bib

@@ -210,3 +210,14 @@ @article{young2016improved
   year={2016},
   URL={https://personal.lse.ac.uk/YoungA/Improved.pdf}
 }
+
+@article{goodnight1980tests,
+  title={Tests of hypotheses in fixed effects linear models},
+  author={Goodnight, James Howard},
+  journal={Communications in Statistics-Theory and Methods},
+  volume={9},
+  number={2},
+  pages={167--180},
+  year={1980},
+  publisher={Taylor \& Francis}
+}

inst/WORDLIST

@@ -1,3 +1,4 @@
+ARMCD
 AstraZeneca
 BCVA
 Benchmarking