diff --git a/NEWS.md b/NEWS.md index 8a1e20190..21432bbbb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,8 @@ ### Bug Fixes - Previously `emmeans` will return `NA` for spatial covariance structure. This is fixed now. +- Previously `car::Anova` will give incorrect results if an interaction term is included and the order of the covariate of interest is not the first categorical variable. This is fixed now. +- Previously `car::Anova` will fail if the model does not contain intercept. 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 diff --git a/R/interop-car.R b/R/interop-car.R index beaafe40f..a007b8bbf 100644 --- a/R/interop-car.R +++ b/R/interop-car.R @@ -40,35 +40,51 @@ h_get_contrast <- function(object, effect, type = c("II", "III", "2", "3"), tol assert_subset(effect, colnames(fcts)) idx <- which(effect == colnames(fcts)) cols <- which(asg == idx) - data <- model.frame(object) - var_numeric <- vapply(data, is.numeric, FUN.VALUE = TRUE) - coef_rows <- length(cols) + xlev <- component(object, "xlev") + contains_intercept <- (!0 %in% asg) && h_first_contain_categorical(effect, fcts, names(xlev)) + coef_rows <- length(cols) - as.integer(contains_intercept) 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)) + if (contains_intercept) { + l_mx[, cols] <- cbind(-1, diag(rep(1, coef_rows))) + } else { + 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]) + additional_numeric <- any(!additional_vars %in% names(xlev)) current_col <- which(asg == i) if (ods[i] >= ods[idx] && all(fcts[, i] >= fcts[, idx]) && !additional_numeric) { sub_mat <- switch(type, "2" = , "II" = { + x1 <- mx[, cols, drop = FALSE] 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 + ret <- solve(t(x1) %*% m %*% x1) %*% t(x1) %*% m %*% x2 + if (contains_intercept) { + ret[-1, ] - ret[1, ] + } else { + ret + } }, "3" = , "III" = { - additional_levels <- vapply(data[additional_vars], function(x) { - if (is.factor(x)) nlevels(x) else length(unique(x)) - }, FUN.VALUE = 1L) - t_levels <- prod(additional_levels) - l_mx[, cols] / t_levels + lvls <- h_obtain_lvls(effect, additional_vars, xlev) + t_levels <- lvls$total + nms_base <- colnames(mx)[cols] + nms <- colnames(mx)[current_col] + nms_base_split <- strsplit(nms_base, ":") + nms_split <- strsplit(nms, ":") + base_idx <- h_get_index(nms_split, nms_base_split) + mt <- l_mx[, cols, drop = FALSE] / t_levels + ret <- mt[, base_idx, drop = FALSE] + # if there is extra levels, replace it with -1/t_levels + ret[is.na(ret)] <- -1 / t_levels + ret } ) l_mx[, current_col] <- sub_mat @@ -114,3 +130,81 @@ Anova.mmrm <- function(mod, type = c("II", "III", "2", "3"), tol = sqrt(.Machine ) ret_df } + + +#' Obtain Levels Prior and Posterior +#' @param var (`string`) name of the effect. +#' @param additional_vars (`character`) names of additional variables. +#' @param xlev (`list`) named list of character levels. +#' @param factors (`matrix`) the factor matrix. +#' @keywords internal +h_obtain_lvls <- function(var, additional_vars, xlev, factors) { + assert_string(var) + assert_character(additional_vars) + assert_list(xlev, types = "character") + nms <- names(xlev) + assert_subset(additional_vars, nms) + if (var %in% nms) { + prior_vars <- intersect(nms[seq_len(match(var, nms) - 1)], additional_vars) + prior_lvls <- vapply(xlev[prior_vars], length, FUN.VALUE = 1L) + post_vars <- intersect(nms[seq(match(var, nms) + 1, length(nms))], additional_vars) + post_lvls <- vapply(xlev[post_vars], length, FUN.VALUE = 1L) + total_lvls <- prod(prior_lvls) * prod(post_lvls) + } else { + prior_lvls <- vapply(xlev[additional_vars], length, FUN.VALUE = 1L) + post_lvls <- 2L + total_lvls <- prod(prior_lvls) + } + list( + prior = prior_lvls, + post = post_lvls, + total = total_lvls + ) +} + +#' Check if the Effect is the First Categorical Effect +#' @param effect (`string`) name of the effect. +#' @param categorical (`character`) names of the categorical values. +#' @param factors (`matrix`) the factor matrix. +#' @keywords internal +h_first_contain_categorical <- function(effect, factors, categorical) { + assert_string(effect) + assert_matrix(factors) + assert_character(categorical) + mt <- match(effect, colnames(factors)) + varnms <- row.names(factors) + # if the effect is not categorical in any value, return FALSE + if (!any(varnms[factors[, mt] > 0] %in% categorical)) { + return(FALSE) + } + # keep only categorical rows that is in front of the current factor + factors <- factors[row.names(factors) %in% categorical, seq_len(mt - 1L), drop = FALSE] + # if previous cols are all numerical, return TRUE + if (ncol(factors) < 1L) { + return(TRUE) + } + col_ind <- apply(factors, 2, prod) + # if any of the previous cols are categorical, return FALSE + return(!any(col_ind > 0)) +} + +#' Test if the First Vector is Subset of the Second Vector +#' @param x (`vector`) the first list. +#' @param y (`vector`) the second list. +#' @keywords internal +h_get_index <- function(x, y) { + assert_list(x) + assert_list(y) + vapply( + x, + \(i) { + r <- vapply(y, \(j) test_subset(j, i), FUN.VALUE = TRUE) + if (sum(r) == 1L) { + which(r) + } else { + NA_integer_ + } + }, + FUN.VALUE = 1L + ) +} diff --git a/man/h_first_contain_categorical.Rd b/man/h_first_contain_categorical.Rd new file mode 100644 index 000000000..d983242c7 --- /dev/null +++ b/man/h_first_contain_categorical.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/interop-car.R +\name{h_first_contain_categorical} +\alias{h_first_contain_categorical} +\title{Check if the Effect is the First Categorical Effect} +\usage{ +h_first_contain_categorical(effect, factors, categorical) +} +\arguments{ +\item{effect}{(\code{string}) name of the effect.} + +\item{factors}{(\code{matrix}) the factor matrix.} + +\item{categorical}{(\code{character}) names of the categorical values.} +} +\description{ +Check if the Effect is the First Categorical Effect +} +\keyword{internal} diff --git a/man/h_get_index.Rd b/man/h_get_index.Rd new file mode 100644 index 000000000..66141806c --- /dev/null +++ b/man/h_get_index.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/interop-car.R +\name{h_get_index} +\alias{h_get_index} +\title{Test if the First Vector is Subset of the Second Vector} +\usage{ +h_get_index(x, y) +} +\arguments{ +\item{x}{(\code{vector}) the first list.} + +\item{y}{(\code{vector}) the second list.} +} +\description{ +Test if the First Vector is Subset of the Second Vector +} +\keyword{internal} diff --git a/man/h_obtain_lvls.Rd b/man/h_obtain_lvls.Rd new file mode 100644 index 000000000..b6d6fa5ff --- /dev/null +++ b/man/h_obtain_lvls.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/interop-car.R +\name{h_obtain_lvls} +\alias{h_obtain_lvls} +\title{Obtain Levels Prior and Posterior} +\usage{ +h_obtain_lvls(var, additional_vars, xlev, factors) +} +\arguments{ +\item{var}{(\code{string}) name of the effect.} + +\item{additional_vars}{(\code{character}) names of additional variables.} + +\item{xlev}{(\code{list}) named list of character levels.} + +\item{factors}{(\code{matrix}) the factor matrix.} +} +\description{ +Obtain Levels Prior and Posterior +} +\keyword{internal} diff --git a/tests/testthat/_snaps/car.md b/tests/testthat/_snaps/car.md index e208a5567..abd9b92a8 100644 --- a/tests/testthat/_snaps/car.md +++ b/tests/testthat/_snaps/car.md @@ -4,7 +4,7 @@ 166.132390274818, 168.517347862691, 148.115819022518, 147.9145531646 ), "F Statistic" = c(36.9114313439307, 0.375714493693334, 31.6630556880233, 142.112225163932, 0.258059683823067), "Pr(>=F)" = c(5.54457484361028e-14, - 0.54074393045656, 7.50169682382158e-08, 2.14434380222247e-43, + 0.54074393045656, 7.5016968238216e-08, 2.14434380222259e-43, 0.855493665048709)), class = c("anova", "data.frame"), row.names = c("RACE", "SEX", "ARMCD", "AVISIT", "ARMCD:AVISIT"), heading = "Analysis of Fixed Effect Table (Type III F tests)") @@ -14,7 +14,7 @@ 174.508554247636, 174.459299901183, 143.379480366266, 143.114582736984 ), "F Statistic" = c(43.8669479408289, 1.16753326235884, 28.5214961211129, 146.699440440163, 0.148199798372577), "Pr(>=F)" = c(3.70573055665334e-16, - 0.281399695369735, 2.86006523470328e-07, 1.67727491048864e-43, + 0.281399695369735, 2.86006523470328e-07, 1.67727491048869e-43, 0.930700208603909)), class = c("anova", "data.frame"), row.names = c("RACE", "SEX", "ARMCD", "AVISIT", "ARMCD:AVISIT"), heading = "Analysis of Fixed Effect Table (Type III F tests)") @@ -24,7 +24,7 @@ 166.132390274818, NA, 168.517347862691, 148.115819022518, 147.9145531646 ), "F Statistic" = c(36.9114313439307, 0.375714493693334, NA, 31.6630556880233, 142.112225163932, 0.258059683823067), "Pr(>=F)" = c(5.54457484361028e-14, - 0.54074393045656, NA, 7.50169682382158e-08, 2.14434380222247e-43, + 0.54074393045656, NA, 7.5016968238216e-08, 2.14434380222259e-43, 0.855493665048709)), class = c("anova", "data.frame"), row.names = c("RACE", "SEX", "SEX2", "ARMCD", "AVISIT", "ARMCD:AVISIT"), heading = "Analysis of Fixed Effect Table (Type III F tests)") diff --git a/tests/testthat/test-car.R b/tests/testthat/test-car.R index 39629e628..a3daa0fa9 100644 --- a/tests/testthat/test-car.R +++ b/tests/testthat/test-car.R @@ -12,6 +12,57 @@ test_that("car emits a message about mmrm registration on load", { expect_message(library(car), "mmrm") }) +# h_first_contain_categorical ---- + +test_that("h_first_contain_categorical works as expected", { + effect <- "FEV1_BL:AVISIT" + factors <- matrix(c(2, 2), nrow = 2L, dimnames = list(c("FEV1_BL", "AVISIT"), c("FEV1_BL:AVISIT"))) + categorical <- c("AVISIT") + expect_true(h_first_contain_categorical(effect, factors, categorical)) + factors2 <- matrix( + c(1, 0, 0, 0, 2, 2), + nrow = 3L, + dimnames = list(c("AGE", "FEV1_BL", "AVISIT"), c("AGE", "FEV1_BL:AVISIT")) + ) + expect_true(h_first_contain_categorical(effect, factors2, categorical)) + categorical2 <- c("NONE") + expect_false(h_first_contain_categorical(effect, factors, categorical2)) + categorical3 <- c("AGE") + expect_false(h_first_contain_categorical(effect, factors, categorical3)) +}) + +# h_obtain_lvls ---- +test_that("h_obtain_lvls works as expected", { + expect_identical( + h_obtain_lvls("a", c("b", "c"), list(a = letters[1:2], b = letters[1:3], c = letters[1:5])), + list( + prior = setNames(integer(0), character(0)), + post = c(b = 3L, c = 5L), + total = 15 + ) + ) + expect_identical( + h_obtain_lvls("a", c("b", "c"), list(b = letters[1:3], c = letters[1:5])), + list( + prior = c(b = 3L, c = 5L), + post = 2L, + total = 15 + ) + ) +}) + +# h_get_index ---- + +test_that("h_get_index works as expected", { + expect_identical( + expect_silent(h_get_index(list(c(1, 2), c(3, 4)), list(c(1), c(3, 4)))), + c(1L, 2L) + ) + expect_identical( + expect_silent(h_get_index(list(c(1, 2), c(1, 4)), list(c(1), c(1, 4)))), + c(1L, NA) + ) +}) # h_get_contrast ---- test_that("h_get_contrast works as expected", { @@ -27,6 +78,140 @@ 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( + 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)) + for (i in seq_len(nrow(ctr))) { + expect_identical( + names(ctr[i, ctr[i, ] != 0]), + sprintf(c("AVISITVIS%s", "RACEBlack or African American:AVISITVIS%s", "RACEWhite:AVISITVIS%s"), i + 1) + ) + } + + 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)) + + for (i in seq_len(nrow(ctr))) { + expect_identical( + names(ctr[i, ctr[i, ] != 0]), + sprintf(c("AVISITVIS%s", "AVISITVIS%s:RACEBlack or African American", "AVISITVIS%s:RACEWhite"), i + 1) + ) + } +}) + +test_that("h_get_contrast works even if only interaction term exists", { + mod1 <- mmrm( + formula = FEV1 ~ FEV1_BL:AVISIT - 1 + ar1(AVISIT | USUBJID), + data = fev_data + ) + ctr <- expect_silent(h_get_contrast(mod1, "FEV1_BL:AVISIT", "3")) + colnames(ctr) <- names(coef(mod1)) + for (i in seq_len(nrow(ctr))) { + expect_identical( + names(ctr[i, ctr[i, ] != 0]), + sprintf(c("FEV1_BL:AVISITVIS%s"), c("1", i + 1)) + ) + } + + mod2 <- mmrm( + formula = FEV1 ~ AVISIT + AVISIT:RACE + FEV1_BL + us(AVISIT | USUBJID), + data = fev_data + ) + ctr2 <- expect_silent(h_get_contrast(mod2, "AVISIT", "3")) + colnames(ctr2) <- names(coef(mod2)) + + for (i in seq_len(nrow(ctr2))) { + expect_identical( + names(ctr2[i, ctr2[i, ] != 0]), + sprintf( + c( + "AVISITVIS%s", "AVISITVIS1:RACEBlack or African American", "AVISITVIS%s:RACEBlack or African American", + "AVISITVIS1:RACEWhite", "AVISITVIS%s:RACEWhite" + ), + i + 1 + ) + ) + } + + mod3 <- mmrm( + formula = FEV1 ~ AVISIT + AVISIT:RACE + FEV1_BL - 1 + us(AVISIT | USUBJID), + data = fev_data + ) + ctr3 <- expect_silent(h_get_contrast(mod3, "AVISIT", "3")) + colnames(ctr3) <- names(coef(mod3)) + + for (i in seq_len(nrow(ctr3))) { + expect_identical( + names(ctr3[i, ctr3[i, ] != 0]), + sprintf( + c( + "AVISITVIS1", "AVISITVIS%s", "AVISITVIS1:RACEBlack or African American", + "AVISITVIS%s:RACEBlack or African American", + "AVISITVIS1:RACEWhite", "AVISITVIS%s:RACEWhite" + ), + as.character(i + 1) + ) + ) + } +}) + +test_that("h_get_contrast works if intercept is not given", { + fit <- mmrm(FEV1 ~ AVISIT * ARMCD - 1 + ar1(AVISIT | USUBJID), data = fev_data) + h_get_contrast(fit, "ARMCD", "2") + h_get_contrast(fit, "AVISIT", "2") + h_get_contrast(fit, "AVISIT:ARMCD", "2") +}) + +test_that("h_get_contrast works for higher-order interaction", { + mod <- mmrm( + formula = FEV1 ~ ARMCD + RACE + AVISIT + RACE * AVISIT * ARMCD + FEV1_BL + ar1(AVISIT | USUBJID), + data = fev_data + ) + ctr <- expect_silent(h_get_contrast(mod, "ARMCD:AVISIT", "3")) + colnames(ctr) <- names(coef(mod)) + for (i in seq_len(nrow(ctr))) { + expect_identical( + names(ctr[i, ctr[i, ] != 0]), + sprintf( + c( + "ARMCDTRT:AVISITVIS%s", + "ARMCDTRT:RACEBlack or African American:AVISITVIS%s", + "ARMCDTRT:RACEWhite:AVISITVIS%s" + ), + as.character(i + 1) + ) + ) + } + + ctr <- expect_silent(h_get_contrast(mod, "RACE:AVISIT", "3")) + colnames(ctr) <- names(coef(mod)) + for (i in seq_len(nrow(ctr))) { + expect_identical( + names(ctr[i, ctr[i, ] != 0]), + sprintf( + c( + "RACE%s:AVISITVIS%s", + "ARMCDTRT:RACE%s:AVISITVIS%s" + ), + c("Black or African American", "White")[(i - 1) %% 2 + 1], + as.character(ceiling(i / 2) + 1) + ) + ) + } }) # Anova ----