fix anova type 3

NEWS.md

@@ -3,6 +3,7 @@
 ### Bug Fixes
 
 - Previously `emmeans` will return `NA` for spatial covariance structure. This is fixed now.
+- Previously `car::Anova` will give incorrect result if the interaction term is included and the order of the covariate of interest is not the first categorical variable. This is fixed now.
 - Previously, `mmrm` will ignore contrasts defined for covariates in the input data set. This is fixed now.
 
 # mmrm 0.3.12

R/interop-car.R

@@ -67,8 +67,17 @@ h_get_contrast <- function(object, effect, type = c("II", "III", "2", "3"), tol
           additional_levels <- vapply(data[additional_vars], function(x) {
             if (is.factor(x)) nlevels(x) else length(unique(x))
           }, FUN.VALUE = 1L)
+          prior_vars <- additional_vars[which(match(additional_vars, row.names(fcts)) < match(effect, row.names(fcts)))]
+          prior_lvls <- vapply(data[prior_vars], function(x) {
+            if (is.factor(x)) nlevels(x) else length(unique(x))
+          }, FUN.VALUE = 1L)
+          prior_lvls <- prod(prior_lvls - 1L)
           t_levels <- prod(additional_levels)
-          l_mx[, cols] / t_levels
+          current_lvls <- length(cols)
+          current_row_idx <- rep(rep(seq_len(current_lvls), each = prior_lvls), times = length(current_col) / prior_lvls / current_lvls)
+          mt <- matrix(0, nrow = current_lvls, ncol = length(current_col))
+          mt[cbind(current_row_idx, seq_len(length(current_row_idx)))] <- 1 / t_levels
+          mt
         }
       )
       l_mx[, current_col] <- sub_mat

tests/testthat/test-car.R

@@ -27,6 +27,50 @@ test_that("h_get_contrast works as expected", {
     h_get_contrast(get_mmrm_trans(), "ARMCD:AVISIT", "3"),
     matrix(rep(rep(c(0, 1), 3), c(6, 1, 9, 1, 9, 1)), nrow = 3, byrow = TRUE)
   )
+  expect_identical(
+    fit <- 
+    h_get_contrast(get_mmrm_trans(), "ARMCD:AVISIT", "3"),
+    matrix(rep(rep(c(0, 1), 3), c(6, 1, 9, 1, 9, 1)), nrow = 3, byrow = TRUE)
+  )
+})
+
+test_that("h_get_contrast works even if the interaction term order changes", {
+  mod1 <- mmrm(
+    formula = FEV1 ~  RACE + AVISIT + RACE * AVISIT + FEV1_BL + us(AVISIT|USUBJID),
+    data = fev_data
+  )
+  ctr <- expect_silent(h_get_contrast(mod1, "AVISIT", "3"))
+  colnames(ctr) <- names(coef(mod1))
+  expect_identical(
+    names(ctr[1, ctr[1, ] != 0]),
+    sprintf(c("AVISITVIS%s","RACEBlack or African American:AVISITVIS%s","RACEWhite:AVISITVIS%s"), "2")
+  )
+  expect_identical(
+    names(ctr[2, ctr[3, ] != 0]),
+    sprintf(c("AVISITVIS%s","RACEBlack or African American:AVISITVIS%s","RACEWhite:AVISITVIS%s"), "3")
+  )
+  expect_identical(
+    names(ctr[3, ctr[3, ] != 0]),
+    sprintf(c("AVISITVIS%s","RACEBlack or African American:AVISITVIS%s","RACEWhite:AVISITVIS%s"), "4")
+  )
+  mod2 <- mmrm(
+    formula = FEV1 ~  AVISIT + RACE + AVISIT * RACE + FEV1_BL + us(AVISIT|USUBJID),
+    data = fev_data
+  )
+  ctr <- expect_silent(h_get_contrast(mod2, "AVISIT", "3"))
+  colnames(ctr) <- names(coef(mod2))
+  expect_identical(
+    names(ctr[1, ctr[1, ] != 0]),
+    sprintf(c("AVISITVIS%s","AVISITVIS%s:RACEBlack or African American","AVISITVIS%s:RACEWhite"), "2")
+  )
+  expect_identical(
+    names(ctr[2, ctr[2, ] != 0]),
+    sprintf(c("AVISITVIS%s","AVISITVIS%s:RACEBlack or African American","AVISITVIS%s:RACEWhite"), "3")
+  )
+  expect_identical(
+    names(ctr[3, ctr[3, ] != 0]),
+    sprintf(c("AVISITVIS%s","AVISITVIS%s:RACEBlack or African American","AVISITVIS%s:RACEWhite"), "4")
+  )
 })
 
 # Anova ----