Small fixes

NEWS.md

@@ -1,5 +1,9 @@
 # mmrm 0.3.18.9000
 
+### Bug Fixes
+
+- `emmeans()` now uses the fitted fixed-effect terms to construct its coefficient basis. This fixes incorrect marginal means and uncertainty estimates for models where the visit variable appears only in interactions, such as treatment-by-visit models without a visit main effect.
+
 ### Miscellaneous
 
 - Satterthwaite, empirical and linear Kenward-Roger covariance preparation now skip unused second derivatives. Linear Kenward-Roger also omits the `R` component and its allocation.

R/interop-emmeans.R

@@ -106,6 +106,13 @@ emm_basis.mmrm <- function(
   grid,
   ...
 ) {
+  # Recovered terms include covariance variables for the reference grid.
+  # Use the fitted fixed-effect terms to preserve the coefficient basis.
+  trms <- stats::delete.response(stats::terms(object))
+  # Be careful with the factor levels:
+  # We only need those that are in the fitted terms.
+  xlev_subset <- intersect(names(xlev), rownames(attr(trms, "factors")))
+  xlev <- xlev[xlev_subset]
   model_frame <- stats::model.frame(
     trms,
     grid,

inst/WORDLIST

@@ -16,6 +16,7 @@ FEV
 GLS
 GmbH
 Gower
+Hadamard
 Hoffmann
 Indexible
 Ingelheim
@@ -50,7 +51,6 @@ ar
 ast
 biostatistics
 blockdiag
-boldsymbol
 cF
 cdot
 cdots
@@ -85,6 +85,7 @@ iG
 ih
 ihj
 ij
+ik
 im
 indexible
 infty

tests/testthat/test-emmeans.R

@@ -366,6 +366,39 @@ test_that("emmeans also works when the visit variable is contained only in an in
   )
 })
 
+test_that("emmeans agrees for equivalent models with and without a visit main effect", {
+  skip_if_not_installed("emmeans", minimum_version = "1.6")
+
+  fit1 <- mmrm(
+    FEV1 ~ RACE + SEX + FEV1_BL + AVISIT + ARMCD:AVISIT + FEV1_BL:AVISIT + us(AVISIT | USUBJID),
+    data = fev_data
+  )
+  fit2 <- mmrm(
+    FEV1 ~ RACE + SEX + FEV1_BL + ARMCD:AVISIT + FEV1_BL:AVISIT + us(AVISIT | USUBJID),
+    data = fev_data
+  )
+
+  # The fits represent the same model, but use different coefficient bases.
+  expect_equal(fitted(fit1), fitted(fit2), tolerance = 1e-6)
+  em1 <- emmeans::emmeans(fit1, ~ ARMCD | AVISIT)
+  em2 <- emmeans::emmeans(fit2, ~ ARMCD | AVISIT)
+  result1 <- as.data.frame(em1)
+  result2 <- as.data.frame(em2)
+
+  expect_equal(result1[c("ARMCD", "AVISIT")], result2[c("ARMCD", "AVISIT")])
+  expect_equal(result1$emmean, result2$emmean, tolerance = 1e-6)
+  expect_equal(result1$SE, result2$SE, tolerance = 1e-6)
+  expect_equal(result1$df, result2$df, tolerance = 1e-6)
+
+  # Build the full design matrix, including aliased columns, from the fitted terms.
+  model_mat <- stats::model.matrix(
+    stats::delete.response(stats::terms(fit2)),
+    stats::model.frame(fit2),
+    contrasts.arg = component(fit2, "contrasts")
+  )
+  expect_identical(colnames(em2@linfct), colnames(model_mat))
+})
+
 test_that("emmeans also works when the visit variable is not part of the covariates", {
   skip_if_not_installed("emmeans", minimum_version = "1.6")
 

vignettes/algorithm.Rmd

@@ -71,7 +71,16 @@ S_i = \begin{pmatrix}
 \]
 $G_i \in \mathbb{R}^{m_i \times m_i}$ is a diagonal matrix with fixed,
 strictly positive weights on its diagonal. It is the identity matrix if no
-weights are specified.
+weights are specified. The `weights` argument to `mmrm()` is a vector with
+one value $w_{ij}$ per observation; $G_i = \operatorname{diag}(w_{i1},
+\dotsc, w_{im_i})$ collects the entries for subject $i$. Thus weights act
+on the variance scale: the marginal residual variance for observation $j$
+is $(S_i^\top \Sigma S_i)_{jj}/w_{ij}$, while covariance between observations
+$j$ and $k$ is divided by $\sqrt{w_{ij}w_{ik}}$. Equivalently, residual
+standard deviations are divided by $\sqrt{w_{ij}}$. For example, weights
+$(2, 1)$ for two observations correspond to $G_i = \operatorname{diag}(2, 1)$;
+the first variance is halved relative to its unweighted value. No normalization
+to a sum of one or to a treatment-group sample size is required.
 Note that this follows from the well known property of the multivariate normal
 distribution that linear combinations of the random vector again have a
 multivariate normal distribution with the correspondingly modified mean
@@ -162,9 +171,9 @@ For example, for an unstructured covariance matrix, $\theta$ has $k = B m(m+1)/2
 #### Spatial covariance matrix
 
 A spatial covariance structure can model individual-specific visit times or
-locations. Let $\boldsymbol{c}_{ij}$ be the coordinate vector of observation
+locations. Let $c_{ij}$ be the coordinate vector of observation
 $j = 1, \dotsc, m_i$ for subject $i$, and define the Euclidean distance
-$d_{i,jr} = \|\boldsymbol{c}_{ij} - \boldsymbol{c}_{ir}\|_2$ between
+$d_{i,jr} = \|c_{ij} - c_{ir}\|_2$ between
 observations $j$ and $r$, with $j,r \in \{1, \dotsc, m_i\}$.
 For one covariance group, the unweighted subject covariance matrix has entries
 \[
@@ -195,6 +204,10 @@ weighted least squares estimator $\hat{\beta}$ solving the estimating equation
 \[
 (X^\top \Omega^{-1} X) \hat{\beta} = X^\top \Omega^{-1} Y.
 \]
+Here $\Omega^{-1}$ is the full inverse residual covariance matrix used by
+weighted least squares (sometimes denoted $W$). It is derived from the
+covariance model and the observation-weight vector through the $G_i$ matrices;
+it is not the vector supplied as `weights`.
 Plugging in $\hat{\beta}$ into the likelihood above gives then the value of the
 function we want to maximize with regards to the variance parameters $\theta$.
 Practically this will be done on the negative log scale:

vignettes/kenward.Rmd

@@ -464,10 +464,11 @@ is a symmetric $c\times c$ matrix. Cyclic invariance of trace gives
   A_2 = \sum_{h,j} W_{hj}\langle K_h,K_j\rangle_F,
 \]
 
-where $\langle\cdot,\cdot\rangle_F$ is the Frobenius inner product.
+where $\langle\cdot,\cdot\rangle_F$ is the [Frobenius inner product](https://en.wikipedia.org/wiki/Frobenius_inner_product).
 With $\mathcal{K}$ containing $\operatorname{vec}(K_h)$ as its $h$th
 column, the implementation calculates
-$A_2 = \operatorname{sum}((\mathcal{K}W)\odot\mathcal{K})$.
+$A_2 = \operatorname{sum}((\mathcal{K}W)\odot\mathcal{K})$, 
+where $\cdot \odot \cdot$ is the [Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)).
 This replaces the pairwise products of $p\times p$ matrices by a dense
 contraction with $c^2$ rows. The remaining scalar formulas above are unchanged.
 The full covariance $W$ is retained, including cross-group entries in grouped
@@ -509,7 +510,7 @@ of $\Sigma$ over our parameters, are non-zero matrices. However, if we use the e
 derivatives are zero matrices. However, the differences are usually small.
 If you would like to match SAS results for the unstructured covariance model, you can use the linear Kenward-Roger approximation.
 
-## Implementations in `mmrm`
+## Implementation in `mmrm`
 
 In package `mmrm`, we have implemented Kenward-Roger calculations based on the previous sections.
 For non-spatial covariance structures, the first-order and second-order derivatives are obtained by
@@ -630,7 +631,7 @@ the inverse-logit reparameterization of $\rho$ produces non-zero second-order
 derivatives where SAS's natural $(0, 1)$-scaled parameterization would give
 zero. The likelihood, $\beta$ estimates and KR-Linear adjusted standard errors
 all match SAS to numerical precision; small (typically sub-percent) differences
-remain on the default Kenward-Roger standard error. To reproduce SAS exactly,
-use the linear Kenward-Roger approximation.
+remain on the default Kenward-Roger standard error. The closest results to SAS are again
+obtain when using the linear Kenward-Roger method in `mmrm`.
 
 # References

vignettes/subsections/_intro-customizations.Rmd

@@ -166,7 +166,15 @@ variable must be coded as a factor.
 
 ## Weighting
 
-Users can perform weighted MMRM by specifying a numeric vector `weights` with positive values.
+Users can perform weighted MMRM by specifying one strictly positive numeric
+weight per observation (row of `data`). These are inverse-variance weights:
+holding the fitted covariance structure fixed, an observation with weight 2
+has half the residual variance of an otherwise comparable observation with
+weight 1 (and its residual standard deviation is divided by $\sqrt{2}$).
+The weights are fixed inputs, not weights on the outcome itself. They do not
+need to sum to one or to any treatment-group sample size; the default is 1
+for every observation. The [algorithm vignette](algorithm.html#linear-model)
+shows how the vector changes the subject covariance matrices.
 
 ```{r common-changes-weights}
 fit_wt <- mmrm(
@@ -177,6 +185,12 @@ fit_wt <- mmrm(
 fit_wt
 ```
 
+Here `fev_data$WEIGHT` supplies varying observation weights. For a simple
+two-level illustration, the same call could use
+`weights = ifelse(fev_data$AVISIT == "VIS1", 2, 1)`. At `VIS1`, this halves
+the model-implied residual variance relative to an otherwise identical
+observation with weight 1.
+
 ## Grouped Covariance Structure
 
 Grouped covariance structures are supported by the`mmrm` package.