diff --git a/NEWS.md b/NEWS.md index 41ed5c4c9..be0494a60 100644 --- a/NEWS.md +++ b/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. diff --git a/R/interop-emmeans.R b/R/interop-emmeans.R index ed8465ee5..eda2d9182 100644 --- a/R/interop-emmeans.R +++ b/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, diff --git a/inst/WORDLIST b/inst/WORDLIST index 4344028e9..80ec11113 100644 --- a/inst/WORDLIST +++ b/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 diff --git a/tests/testthat/test-emmeans.R b/tests/testthat/test-emmeans.R index 5e99fc4da..e01770f4b 100644 --- a/tests/testthat/test-emmeans.R +++ b/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") diff --git a/vignettes/algorithm.Rmd b/vignettes/algorithm.Rmd index 99aa1c8e9..6f6d062cd 100644 --- a/vignettes/algorithm.Rmd +++ b/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: diff --git a/vignettes/kenward.Rmd b/vignettes/kenward.Rmd index f29d58573..76e9de0b6 100644 --- a/vignettes/kenward.Rmd +++ b/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 diff --git a/vignettes/subsections/_intro-customizations.Rmd b/vignettes/subsections/_intro-customizations.Rmd index 299907ad3..f5bbb85fa 100644 --- a/vignettes/subsections/_intro-customizations.Rmd +++ b/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.