Add G-computation correction for emmeans, contrasts argument, and expanded emp_start

.Rbuildignore

@@ -21,3 +21,4 @@ simulations
 ^data-raw$
 ^cache$
 ^[.]?air[.]toml$
+^\.claude$

.gitignore

@@ -1,3 +1,4 @@
+.claude
 .DS_Store
 .httr-oauth
 .project

.pre-commit-config.yaml

@@ -3,7 +3,7 @@
 # R specific hooks: https://github.com/lorenzwalthert/precommit
 repos:
   - repo: https://github.com/lorenzwalthert/precommit
-    rev: v0.4.3.9003
+    rev: v0.4.3.9021
     hooks:
       - id: style-files
         args:
@@ -89,7 +89,7 @@ repos:
           README.md
           )$
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v5.0.0
+    rev: v6.0.0
     hooks:
       - id: check-added-large-files
         args: ["--maxkb=200"]

DESCRIPTION

@@ -1,11 +1,12 @@
 Type: Package
 Package: mmrm
 Title: Mixed Models for Repeated Measures
-Version: 0.3.17.9000
+Version: 0.3.17.9005
 Authors@R: c(
     person("Daniel", "Sabanes Bove", , "daniel.sabanes_bove@rconis.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-0176-9239")),
-    person("Liming", "Li", , "liming.li1@astrazeneca.com", role = "aut", comment = c(ORCID = "0009-0008-6870-0878")),
+    person("Liming", "Li", , "liming.li1@astrazeneca.com", role = "aut",
+           comment = c(ORCID = "0009-0008-6870-0878")),
     person("Julia", "Dedic", , "julia.dedic@roche.com", role = "aut"),
     person("Doug", "Kelkhoff", , "doug.kelkhoff@roche.com", role = "aut"),
     person("Kevin", "Kunzmann", , "kevin.kunzmann@boehringer-ingelheim.com", role = "aut"),
@@ -108,6 +109,7 @@ NeedsCompilation: yes
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.3
 Collate:
+    'gcomp.R'
     'between-within.R'
     'catch-routine-registration.R'
     'component.R'

NEWS.md

@@ -1,3 +1,11 @@
+# mmrm 0.3.17.9005
+
+### New Features
+
+- `mmrm()` and `fit_mmrm()` now accept a `contrasts` argument, allowing users to specify contrast matrices or functions for factor variables, matching the interface of `lm()`. When an explicit contrast matrix includes levels not present in the fitting data, those levels are preserved in the model and marked as aliased, enabling prediction on new data containing those levels (#562).
+- `emp_start()` now supports all non-spatial covariance structure types, not just unstructured. Parametrized covariance structures now use a simple moment-matching procedure to identify a reasonable empirical starting value.
+- `mmrm_control()` gains `emmeans_gcomp_vars` argument, enabling G-computation correction in `emmeans()` output for models with covariate-by-treatment interactions. When set, `emmeans()` returns the average treatment effect (ATE) with standard errors that account for covariate variability across subjects. See `?emmeans_support` for details.
+
 # mmrm 0.3.17.9000
 
 ### Bug Fixes

R/fit.R

@@ -222,6 +222,14 @@ refit_multiple_optimizers <- function(fit, ..., control = mmrm_control(...)) {
 #' @param disable_theta_vcov (`flag`)\cr whether to disable calculation of
 #'   variance-covariance matrix for variance parameters. This can speed up fitting
 #'   when there are many variance parameters, see details.
+#' @param emmeans_gcomp_vars (`character` or `NULL`)\cr names of variables to
+#'   treat as fixed in the G-computation correction for emmeans. When non-`NULL`,
+#'   enables the correction in [`emmeans_support`], producing the average
+#'   treatment effect (ATE) with corrected standard errors. The visit variable
+#'   (identified from the covariance structure) is computed separately per
+#'   level. All other named variables are treated as intervention (subjects
+#'   pooled across levels). All variables not named are averaged over using
+#'   each subject's actual values. Defaults to `NULL` (no correction).
 #' @param ... additional arguments passed to [h_get_optimizers()].
 #'
 #' @details
@@ -282,11 +290,13 @@ mmrm_control <- function(
   accept_singular = TRUE,
   drop_visit_levels = TRUE,
   disable_theta_vcov = FALSE,
+  emmeans_gcomp_vars = NULL,
   ...,
   optimizers = h_get_optimizers(...)
 ) {
   assert_count(n_cores, positive = TRUE)
   assert_character(method)
+  assert_character(emmeans_gcomp_vars, min.len = 1L, null.ok = TRUE)
   if (is.null(start)) {
     start <- std_start
   }
@@ -352,7 +362,8 @@ mmrm_control <- function(
       vcov = vcov,
       n_cores = as.integer(n_cores),
       drop_visit_levels = drop_visit_levels,
-      disable_theta_vcov = disable_theta_vcov
+      disable_theta_vcov = disable_theta_vcov,
+      emmeans_gcomp_vars = emmeans_gcomp_vars
     ),
     class = "mmrm_control"
   )
@@ -374,6 +385,16 @@ mmrm_control <- function(
 #'   as produced with [cov_struct()], or value that can be coerced to a
 #'   covariance structure using [as.cov_struct()]. If no value is provided,
 #'   a structure is derived from the provided formula.
+#' @param contrasts (`list` or `NULL`)\cr an optional named list of contrast
+#'   matrices or contrast functions (like [stats::contr.sum] or
+#'   [stats::contr.poly]) for specific factor variables, matching the
+#'   `contrasts` argument in [stats::lm()]. The list names must correspond to
+#'   factor variable names in the model formula. When `NULL` (the default),
+#'   the contrasts set on the factor variables in `data` are used. If a
+#'   contrast matrix has rownames that include levels not present in `data`,
+#'   those levels are preserved and the corresponding model matrix columns
+#'   are marked as aliased (not estimable), enabling prediction on new data
+#'   containing those levels.
 #' @param control (`mmrm_control`)\cr fine-grained fitting specifications list
 #'   created with [mmrm_control()].
 #' @param ... arguments passed to [mmrm_control()].
@@ -466,13 +487,16 @@ mmrm <- function(
   data,
   weights = NULL,
   covariance = NULL,
+  contrasts = NULL,
   reml = TRUE,
   control = mmrm_control(...),
   ...
 ) {
   assert_false(!missing(control) && !missing(...))
   assert_class(control, "mmrm_control")
   assert_list(control$optimizers, min.len = 1)
+  assert_list(contrasts, null.ok = TRUE, names = "unique")
+  emmeans_gcomp_vars <- control$emmeans_gcomp_vars
 
   if (control$method %in% c("Kenward-Roger", "Kenward-Roger-Linear") && !reml) {
     stop("Kenward-Roger only works for REML")
@@ -494,14 +518,21 @@ mmrm <- function(
   } else {
     attr(weights, which = "dataname") <- deparse(match.call()$weights)
   }
+
+  # Validate G-computation fixed variables if specified.
+  if (!is.null(emmeans_gcomp_vars)) {
+    assert_subset(emmeans_gcomp_vars, names(data))
+  }
+
   tmb_data <- h_mmrm_tmb_data(
     formula_parts,
     data,
     weights,
     reml,
     singular = if (control$accept_singular) "drop" else "error",
     drop_visit_levels = control$drop_visit_levels,
-    allow_na_response = FALSE
+    allow_na_response = FALSE,
+    contrasts = contrasts
   )
   fit <- structure("", class = "try-error")
   names_all_optimizers <- names(control$optimizers)
@@ -589,6 +620,20 @@ mmrm <- function(
     stop("Unrecognized coefficent variance-covariance method!")
   }
 
+  # G-computation correction: store metadata for emmeans hook.
+  fit$emmeans_gcomp_vars <- emmeans_gcomp_vars
+  if (!is.null(emmeans_gcomp_vars)) {
+    # Store subject-level covariate data from the ORIGINAL data (pre-NA-removal).
+    # This includes subjects with observed covariates but missing outcomes,
+    # who contribute to the G-computation average but not to beta estimation.
+    subject_var <- formula_parts$subject_var
+    subj_rows <- !duplicated(data[[subject_var]])
+    model_vars <- all.vars(formula_parts$model_formula)
+    keep_vars <- intersect(c(subject_var, model_vars), names(data))
+    fit$emmeans_gcomp_subject_data <- data[subj_rows, keep_vars, drop = FALSE]
+    rownames(fit$emmeans_gcomp_subject_data) <- NULL
+  }
+
   class(fit) <- c("mmrm", class(fit))
   fit
 }