---
title: "Adding anova() to mmrm"
author: "Nik Krieger"
format: gfm
editor: visual
---

## Setup

```{r setup}
library(mmrm)
```

# Background

A current drawback of `mmrm` is that it does not support `stats::anova()` (i.e., `anova.mmrm()` does not exist). This is discussed here: <https://github.com/openpharma/mmrm/issues/164>. Users are interested in an `mmrm` method for `anova()` supporting the likelihood ratio rest (LRT) for `mmrm` model fits.

# Specifications

## Primary Usage: LRT

The primary usage of `anova.mmrm()` would be to compare multiple nested model fits with a focus on performing LRT on adjacent models. To accomplish this, supplied models must meet the following requirements.

### Compatible Covariates

If all supplied models were estimated using maximum likelihood (ML), the models must have nested covariates in order to undergo LRT. In other words, the set of covariates for each model must be a subset of the covariates of the next model. See: <https://github.com/openpharma/mmrm/issues/164#issuecomment-1446364890>.

However, if any of the supplied models were estimated using restricted maximum likelihood (REML), all models must have the same covariates.

### Nested covariance structures

LRT requires that the covariance structure of each model must be either (1) the same as that of the next model or (2) a special case of that of the next model. For instance, the homogeneous first-order autoregressive structure (`ar1`) is a special case of the Toeplitz structure (`toep`), so models using `ar1` must precede models using `toep`. Concerning models with unstructured covariances, they must follow all models having a covariance structure. See: <https://github.com/openpharma/mmrm/issues/164#issuecomment-1446627674> as well as the **Nested covariance structures proofs** appendix below.

### **Increasing degrees of freedom**

Each supplied model fit must have more degrees of freedom (*df*) than the preceding model fit.

### Common data sets

LRT compares models that are fit to the same data. Model fits must either already use the same data or be refitted to a set of common observations in their respective data sets. See <https://github.com/openpharma/mmrm/issues/164#issuecomment-1446839199> and <https://github.com/openpharma/mmrm/issues/164#issuecomment-1447419569>.

## Compliance with `?stats::anova` 

As a method of the the generic function `stats::anova()`, `anova.mmrm()` will comply with the assertions in the generic's help page, `?stats::anova`.

-   This method will yield a table of results.

-   This method will accept a sequence of `mmrm` fits that are tested against one another in the order supplied.

-   When a single `mmrm` fit is supplied, the resulting table will test the significance of the model's terms.

## Following the lead of `anova()` methods in `nlme` and `lme4`

In order that `anova.mmrm()` may be comparable to `anova()` methods in the `nlme` and `lme4` packages, it shall deliver a table of fit diagnostics (e.g., AIC, BIC, log likelihood) regardless of whether LRT is performed.

Furthermore, `anova.mmrm()` shall accept the following `TRUE`/`FALSE` arguments:

-   `test` – opt out of obtaining LRT results (see `?nlme::anova.lme`).

-   `refit` – refit the models if they are incomparable (see `?lme4::anova.merMod`).

# Prerequisite: augment `logLik.mmrm_tmb()`

This function shall return more data in the `attributes()` of the returned value:

-   `n_param`: the number of variance parameters.
-   `n_coef`: the number of estimated coefficients.
-   `df`: the degrees of freedom
    -   If REML estimation was used, it shall equal `n_param`.
    -   If ML estimation was used, it shall equal `n_param + n_coef`.
-   Add `"logLik"` to the `class()` of the returned value

Prototype:

```{r}
#' @describeIn mmrm_tmb_methods obtains the attained log likelihood value.
#' @importFrom stats logLik
#' @exportS3Method
#' @examples
#' # Log likelihood given the estimated parameters:
#' logLik(object)
logLik.mmrm_tmb <- function(object, ...) {
  out <- -component(object, "neg_log_lik")

  # Number of parameters
  n_param <- component(object, "n_theta")

  # Number of estimated coefficients
  n_coef <- length(coef(object, complete = FALSE))
  
  df <- n_param + n_coef * !component(object, "reml")
  
  attr(out, "n_param") <- n_param
  attr(out, "n_coef") <- n_coef
  attr(out, "df") <- df
  class(out) <- "logLik"
  
  out
}
```

# Designing `anova.mmrm()`

## Function specification

### Usage

```{r eval=FALSE}
## S3 method for class 'mmrm'
anova(object, ..., test = TRUE, refit = FALSE)
```

### Arguments

|  |  |
|------------------------------------|------------------------------------|
| `object` | An object inheriting from class `mmrm`: a fitted mixed model with repeated measures. |
| `...` | Other optional fitted model objects inheriting from class `mmrm`. |
| `test` | Logical indicating whether the output should include likelihood ratio test (LRT) results comparing the model fits to one another. Defaults to `TRUE`. |
| `refit` | Logical indicating whether the models should be refitted with the dataset consisting of their shared set of observations before performing diagnostics and testing. This is ignored if the models already share the same dataset. If `refit = FALSE` and the models have different underlying data sets, an error will be thrown. Defaults to `FALSE`. |

## Function outline

### Plain language

#### Special case: supplying only one model fit

If the user provides only one model fit to `anova.mmrm()`, the function will calculate and return the significance of the model terms in a data frame.

#### Comparing multiple model fits

If the user provides more than one model fit, a data frame is produced with a row for each model, and the columns will contain identifying information as well as standard diagnostics (e.g., AIC; log likelihood).

If the user supplied `test = TRUE`, the model fits are subject to further tests to ensure they are nestable. If the datasets for the fits are not the same and the user specified `refit = TRUE`, the observations common to all the fits' datasets are compiled, and all models are refitted using this data set. The pairs of models then undergo likelihood ratio testing, and the results are augmented into the data frame of results.

### Prototype code

#### `anova.mmrm()`

```{r }
anova.mmrm <- function (object, ..., test = TRUE, refit = FALSE) {
  
  checkmate::assert_class(object, "mmrm")
  
  fits <- list(object, ...)
  # If there are objects supplied to "..."
  if (length(fits) == 1L) {
    out <- h_anova_single_mmrm_model(object)
  } else {
    
    # Ensure all objects in ... are mmrm fits
    lapply(fits[-1L], checkmate::assert_class, classes = "mmrm", .var.name = "...")
    
    # If the data are not all the same and refit = TRUE, refit all with the 
    # largest common data set
    if (refit && !h_check_fits_all_data_same(fits)) {
      common_data <- h_fits_common_data(fits)
      needs_refit <- 
        !vapply(
          lapply(fits, h_get_minimal_fit_data), 
          h_check_columns_nested,
          FUN.VALUE = logical(1L),
          data_augmented = common_data
        )
      fits[needs_refit] <-
        lapply(fits[needs_refit], h_refit_mmrm, data = common_data)
    } else {
      needs_refit <- rep_len(FALSE, length(fits))
    }
    
    # vector of log likelihood results for the fits
    logLik_vec <- lapply(fits, logLik.mmrm_tmb)
    
    # Calculate the standard diagnostics.
    out <-
      data.frame(
        Model = seq_along(fits),
        refit = needs_refit,
        REML = vapply(fits, component, logical(1L), "reml"),
        n_param = vapply(logLik_vec, attr, numeric(1L), "n_param"),
        n_coef = vapply(logLik_vec, attr, numeric(1L), "n_coef"),
        df = vapply(logLik_vec, attr, numeric(1L), "df"), 
        AIC = vapply(fits, AIC, numeric(1)),
        BIC = vapply(fits, BIC, numeric(1)),
        logLik = as.numeric(logLik_vec)
      )
    
    # If the user has requested the likelihood ratio test...
    if (test) {
      
      h_assert_LRT_suitability(fits, refit, dfs = out$df, REMLs = out$REML)
        
      # Label for the pair being compared (e.g., "3 vs 4")
      out$test <- c(NA, paste(out$Model[-length(fits)], "vs", out$Model[-1L]))
      
      # Calculate the log of likelihood ratios between the pairs using
      # base::diff()
      out$log_likelihood_ratio <- c(NA, diff(out$logLik))
      
      # Calculate the p-value using the log likelihood ratios and differences in
      # df
      out$p_value <-
        pchisq(
          2 * abs(out$log_likelihood_ratio),
          df = abs(diff(out$df)), 
          lower.tail = FALSE
        )
    }
    
    out$call <- vapply(lapply(fits, component, "call"), deparse1, character(1L))
    
  }
  
  class(out) <- union("anova.mmrm", class(out))
  
  out
}
```

#### Supporting functions:

```{r}
# Checks that the underlying data sets for all the fits are the same. Before
# comparing, it assembles the minimal set of columns actually used for each fit.
# The following three criteria must be met for this function to return TRUE:
# 1) All data sets must have the same number of rows.
# 2) Each model's covariates must be a subset of the next model's covariates.
# 3) The columns in common between neighboring data sets must be identical.
h_check_fits_all_data_same <- function(fits) {
  # Grab the data frames behind the model fits into a separate list.
  datasets <- lapply(fits, h_get_minimal_fit_data)
  
  # Ensure the number of rows of all data sets are the same.
  if (length(unique(vapply(datasets, nrow, numeric(1L)))) > 1L) {
    return(FALSE)
  }
  
  cols <- lapply(datasets, colnames)
  
  for (i in seq_along(fits)[-1L]) {
    # Ensure previous dataset's cols are a subset of the current dataset's cols
    if (anyNA(match(cols[[i - 1L]], cols[[i]]))) {
      return(FALSE)
    }
  }
  
  # Sort the first data frame.
  datasets[[1L]] <- h_dataset_sort_all(datasets[[1L]])
  for (i in seq_along(fits)[-1L]) {
    
    # Pull prev dataset's columns to the front of current dataset. Then sort.
    cols[[i]] <- union(cols[[i - 1L]], cols[[i]])
    datasets[[i]] <- h_dataset_sort_all(datasets[[i]][cols[[i]]])
    
    # Ensure the common columns are the same
    if (
      !isTRUE(all.equal(
        datasets[[i - 1L]], 
        datasets[[i]][cols[[i - 1L]]],
        check.attributes = FALSE
      ))
    ) {
      return(FALSE)
    }
  }
  
  TRUE
}



# Iteratively merge the data sets underlying the fits to get the common rows and
# all columns.
h_fits_common_data <- function(fits) {
  datasets <- lapply(fits, h_get_minimal_fit_data)
  out <- Reduce(merge, datasets)
  out <- droplevels(out)
  out
}



# Predicate function indicating whether one dataset is a subset of a second
# dataset's columns. Returns TRUE if the subset of the second dataset's columns
# results in the first dataset.
h_check_columns_nested <- function(data_basic, data_augmented) {
  if (nrow(data_basic) != nrow(data_augmented)) {
    return(FALSE)
  }
  
  colnames_basic <- colnames(data_basic)
  
  if (anyNA(match(colnames_basic, colnames(data_augmented)))) {
    return(FALSE)
  }
  
  isTRUE(all.equal(
    h_dataset_sort_all(data_basic),
    h_dataset_sort_all(data_augmented[colnames_basic]),
    check.attributes = FALSE
  ))
}




# Refit a single mmrm model with a new dataset
h_refit_mmrm <- function(fit, data) {
  
  # Grab the model call
  expr <- component(fit, "call")
  
  # Grab the environment of the model's formula
  env <- environment(fit[["formula_parts"]][["full_formula"]])
  
  # Generate a name guaranteed not to be in env nor its enclosing frames
  data_name <- h_generate_new_name("data", env)
  
  # Create a child environment whose parent is the model's environment
  env <- new.env(parent = env)
  
  # bind the new data to the new name in the new, child environment
  env[[data_name]] <- data
  
  # Force the data argument of the model call to be the new name
  expr[["data"]] <- as.name(data_name)
  
  # Evaluate the updated model call in the appropriate environment
  fit <- eval(expr, env)
  
  fit
}



# Generates a name that is guaranteed not to be found in env nor its enclosing
# frames
h_generate_new_name <- function(x, env) {
  # first, force x to be a syntactically valid name(s)
  x <- make.names(x, unique = TRUE)
  
  # As long as we keep finding an object in env (or its enclosing frames) whose
  # name is the same as the last element of x...
  while (exists(x[length(x)], envir = env, inherits = TRUE)) {
    # ...copy the first element of x and append it to the end. Run
    # make.names(unique = TRUE) again in the hopes that the resulting last
    # element will not be found on the search path next time
    x <- make.names(x[c(seq_along(x), 1L)], unique = TRUE)
  }
  # If we've gotten here, the last element of doesn't exist in env or its
  # enclosing frames
  x[length(x)]
}



# Function that checks the fits to make sure LRT is appropriate. Checks degrees
# of freedom, nesting of covariates (in light of ML/REML usage), and nesting of
# covariance structures.
h_assert_LRT_suitability <- function(fits, refit, dfs, REMLs) {
  
  if (any(diff(dfs) <= 0)) {
    stop(
      "The degrees of freedom (df) of each model must be less than the df of ",
      "the next model in order to perform likelihood ratio testing (LRT). ",
      "Bypass LRT with test = FALSE.",
      call. = FALSE
    )
  }
  
  # Do any models use REML?
  any_reml <- any(REMLs)
  
  # First iterate through and check the models to make sure their covariates
  # and covariance structures are nested
  for (i in seq_along(fits)[-1L]) {
    h_assert_nested_models(fits[[i - 1L]], fits[[i]], any_reml = any_reml)
  }
  
  # If we didn't refit, throw errors if not all data are the same
  if (!refit && !h_check_fits_all_data_same(fits)) {
    stop(
      "Likelihood ratio testing requires all fits to use the same data.",
      " Not all fits have the same data and refit = FALSE.",
      " Consider setting test = FALSE or refit = TRUE.",
      call. = FALSE
    )
  }
  
  TRUE
}



# Compares two models to see if they are properly nested.
# Errors out if they don't have the same visit, subject, and grouping vars
# Errors out if the first model's covariates isn't a subset of the second
# Errors out if the first model's covariance structure isn't a special case of
# the second model's covaraince structure.
# Warns if the models have the same covariates and covariance structure
h_assert_nested_models <- function(model_basic, model_augmented, any_reml) {
  
  model_basic <- model_basic[["formula_parts"]]
  model_augmented <- model_augmented[["formula_parts"]]

  if (!identical(model_basic[["visit_var"]], model_augmented[["visit_var"]])) {
    stop("Models must all have the same visit variable.", call. = FALSE)
  }

  if (!identical(model_basic[["subject_var"]], model_augmented[["subject_var"]])) {
    stop("All models must have the same subject variable.", call. = FALSE)
  }

  if (!identical(model_basic[["group_var"]], model_augmented[["group_var"]])) {
    stop("All models must have the same grouping variable.", call. = FALSE)
  }
  
  covar_nesting <- h_check_covar_nesting(model_basic, model_augmented)
  if (any_reml && covar_nesting == "nested") {
    stop("If any models use REML, all models' covariates must be the same.",
         call. = FALSE)
  }
  
  cov_struct_nesting <- h_check_cov_struct_nesting(model_basic, model_augmented)
  
  if (covar_nesting == "identical" && cov_struct_nesting == "identical") {
    warning("Two models in the sequence have identical covariates and ",
            "covariance structures.", call. = FALSE)
  }
  
  NULL
}



# Checks the covariates of two models to see if they are nested (or identical).
# The first model's covariates must be a subset of the second model's covariates
h_check_covar_nesting <- function(model_basic, model_augmented) {
  basic_covars <- attr(terms(model_basic[["model_formula"]]), "term.labels")
  aug_covars <- attr(terms(model_augmented[["model_formula"]]), "term.labels")

  if (anyNA(match(basic_covars, aug_covars))) {
    stop("Each model's covariates must be a subset of the next model's ",
         "covariates.", call. = FALSE)
  }

  if (anyNA(match(aug_covars, basic_covars))) {
    "nested"
  } else {
    "identical"
  }
}



# Checks the covariance structures of two models to see if they are nested (or
# identical). This is a check against a mathematically determined hierarchy.
h_check_cov_struct_nesting <- function(model_basic, model_augmented) {
  basic_cov_struct <- model_basic[["cov_type"]]
  aug_cov_struct <- model_augmented[["cov_type"]]
  
  # Hierarchy, specifying the "parent" covariance structures of each available 
  # type in mmrm
  cov_struct_nesting <-
    list(
      ad = c("adh", "us"),
      adh = "us",
      ar1 = c("ad", "adh", "ar1h", "toep", "toeph", "us"),
      ar1h = c("adh", "toeph", "us"),
      cs = c("csh", "toep", "toeph", "us"),
      csh = c("toeph", "us"),
      toep = c("toeph", "us"),
      toeph = "us",
      us = character(0L),
      sp_exp = character(0L)
    )
  
  if (identical(basic_cov_struct, aug_cov_struct)) {
    "identical"
  } else if (any(cov_struct_nesting[[basic_cov_struct]] == aug_cov_struct)) {
    "nested"
  } else {
    stop("Each model's covariance structure must be either identical to or a ",
         "special case of the next model's covariance structure.",
         call. = FALSE)
  }
}



# Returns a character vector of the names of the columns actually used in a fit
h_get_minimal_fit_data <- function(fit) {
  predictors <- 
    unique(
      unlist(
        fit[["formula_parts"]][
          c("visit_var", "subject_var", "group_var", "model_var")
        ],
        use.names = FALSE
      )
    )
  terms_attr <- attributes(terms(fit))
  response <- as.character(terms_attr$variables[[terms_attr$response + 1]])
  fit[["data"]][c(response, predictors)]
}



# Sorts a data frame by its columns' values in the order the columns appear
h_dataset_sort_all <- function(data) {
  data[do.call(order, unname(data)), ]
}


# A function that creates the table of `anova.mmrm()` results when only one
# model object is provided. This table contains the significance of each of the
# model's terms.
h_anova_single_mmrm_model <- function(object) {
  # attributes of terms.object for the model formula without the cov struct
  terms_attr <- attributes(terms(object$formula_parts$model_formula))
  
  # Get a character vector of the terms.
  # Prepend with "(Intercept)" if an intercept is present
  terms <- c(if (terms_attr$intercept) "(Intercept)", terms_attr$term.labels)
  
  # Vector identifying each coefficient's corresponding term
  assign_vec <- attr(component(object, "x_matrix"), "assign")
  
  # Create an identity matrix with a row/col for each coefficient
  identity_mtx <- diag(length(assign_vec))
  
  # Split the rows of identity_mtx according to assign_vec, creating a list of
  # contrast matrices: one for each term in the model.
  contrasts <-
    split.data.frame(identity_mtx, factor(assign_vec, labels = terms))
  
  # Get df_md() results for each contrast matrix. Result is a list of lists
  df_md_results <- lapply(contrasts, df_md, object = object)
  
  # Turn results into a data frame. The row.names identify the terms.
  out <- do.call(rbind, lapply(df_md_results, as.data.frame))
  
  out
}
```

# Examples

## Setup

Using `mmrm::fev_data`, here are three models: one including `SEX`, one including `SEX` and `RACE`, and one including `SEX` and `RACE` and the interaction between the two. All three use the `ar1` covariance structure, and all three use maximum likelihood (ML) estimation, as specified by `reml = FALSE`.

```{r}
fit_sex <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX + ARMCD + ar1(AVISIT | USUBJID),
    data = fev_data,
    reml = FALSE
  )

fit_sex_race <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX + RACE + ARMCD + ar1(AVISIT | USUBJID),
    data = fev_data,
    reml = FALSE
  )

fit_interaction <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX * RACE + ARMCD + ar1(AVISIT | USUBJID),
    data = fev_data,
    reml = FALSE
  )
```

## Single model

First, run `anova()` on just one model fit:

```{r}
anova(fit_interaction)
```

This yields the significance of each term in the model.

## Multiple models

Now, use multiple models. Note that the likelihood ratio test is included by default.

```{r}
anova(fit_sex, fit_sex_race, fit_interaction)
```

### Nesting checks

#### `df` must increase

Models can only increase in *df* from left to right. Thus, successive models can only gain terms from left to right, so changing the order of the models in the code above will result in an error.

```{r error=TRUE}
anova(fit_sex, fit_interaction, fit_sex_race)
```

As the error message states, we can obtain results without LRT—which still includes *df*—by specifying `test = FALSE`:

```{r}
anova(fit_sex, fit_interaction, fit_sex_race, test = FALSE)
```

Notice that the *df* from the second to the third models, `fit_interaction` and `fit_sex_race`, decreases from 10 to 8, so `anova.mmrm()` will not perform LRT.

#### Nested covariates

ML-estimated models that pass the increasing *df* check have to subsequently pass the nested covariates check: the terms in each model must be a subset of the terms in the next model. Note that if at least one model estimated with restricted maximum likelihood (REML) is present, all covariates in all models must be identical (see the **Greater Restrictions when using REML estimation** section below).

Let's create a new model, `fit_race`, that is like `fit_sex` with `RACE` substituted for `SEX`. Because `RACE` has three levels as opposed to the two levels in `SEX`, it has more *df*. We can see this by adding it the call to `anova(test = FALSE)`:

```{r}
fit_race <- 
  mmrm(
    FEV1 ~ FEV1_BL + RACE + ARMCD + ar1(AVISIT | USUBJID),
    data = fev_data,
    reml = FALSE
  )

anova(fit_sex, fit_race, fit_sex_race, fit_interaction, test = FALSE)
```

The four models in the `anova()` call above have increasing *df*, but LRT still cannot be performed:

```{r error=TRUE}
anova(fit_sex, fit_race, fit_sex_race, fit_interaction)
```

`fit_sex` has fewer *df* than `fit_race` but its covariates are not nested: `fit_race` does not include `SEX`.

#### Nested covariance structures

The covariance structures must also be nested from model to model. Each model's covariance structure must be identical to or a special case of the next model's covariance structure.

Let's create three new models, `fit_sex_csh`, `fit_sex_toeph`, and `fit_sex_us`, which are like `fit_sex` but with the `csh`, `toeph`, and `us` covariance structures instead of `ar1`. Note that `ar1` and `csh` are both special cases of `toeph` and `us`, so `fit_sex` (which uses `ar1`) and `fit_sex_csh` are each considered nested within `fit_sex_toeph` and `fit_sex_us`; however, `fit_sex` nor `fix_sex_csh` are special cases of the other, so they cannot be nested.

```{r}
fit_sex_csh <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX + ARMCD + csh(AVISIT | USUBJID),
    data = fev_data,
    reml = FALSE
  )

fit_sex_toeph <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX + ARMCD + toeph(AVISIT | USUBJID),
    data = fev_data,
    reml = FALSE
  )

fit_sex_us <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX + ARMCD + us(AVISIT | USUBJID),
    data = fev_data,
    reml = FALSE
  )
```

Let's first view the *df* of these models, which all have identical covariates, using `anova(test = FALSE)`:

```{r}
anova(fit_sex, fit_sex_csh, fit_sex_toeph, fit_sex_us, test = FALSE)
```

However, despite their increasing *df* and identical covariates, we cannot turn LRT testing on with `test = TRUE` because `fit_sex` is not nested within `fit_sex_csh`:

```{r error=TRUE}
anova(fit_sex, fit_sex_csh, fit_sex_toeph, fit_sex_us, test = TRUE)
```

If we follow the covariance structure nesting requirements, we can perform LRT with different subsets of these four models:

```{r}
anova(fit_sex, fit_sex_toeph, fit_sex_us, test = TRUE)
anova(fit_sex_csh, fit_sex_toeph, fit_sex_us, test = TRUE)
anova(fit_sex, fit_sex_us, test = TRUE)
anova(fit_sex_csh, fit_sex_us, test = TRUE)
```

### Data sameness checks and the `refit` argument

#### Setup

First, let's refit `fev_sex_toeph` based on a subset of the `mmrm::fev_data` dataset such that each subject only has three time points instead of 4:

```{r}
fev_subset <- droplevels(fev_data[fev_data$VISITN < 4, ])

fit_sex_toeph_subset <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX + ARMCD + toeph(AVISIT | USUBJID),
    data = fev_subset,
    reml = FALSE
  )

fev_subset
```

Without performing LRT, let's first inspect the original `fit_sex` and the new `fit_sex_toeph_subset` along with `fit_sex_toeph` and `fit_sex_us`:

```{r}
anova(fit_sex, fit_sex_toeph_subset, fit_sex_toeph, fit_sex_us, test = FALSE)
```

Even though the `toeph` models (Models 2 and 3) have the same covariates and covariance structures, we see that they have different *df* because the number of time points can affect the number of variance parameters depending on the covariance structure that was used; `toeph` is one such covariance structure. For ML-estimated models in particular:

$$
\text{df} = \text{number of coefficients} + \text{number of variance parameters}
$$

For all these models, there are four coefficients (one for each of the three predictors and one for the intercept). For the `toeph` covariance structure, the number of variance parameters is $2m - 1$, where $m$ is the number of time points (see `?mmrm::cov_types`).

Therefore, we can verify that the original `fit_sex_toeph` model's *df* is $4 + 2\times 4 - 1 = 11$, and the new `fit_sex_toeph_subset` model's *df* is $4 + 2\times 3 - 1 = 9$.

#### Data sameness assertion: Different data, `refit = FALSE`

Because the *df* increases with each model, the covariates are nested, and the covariance structures are also nested, it may appear at first that we can perform the following likelihood ratio testing as is; however, this will not work, because LRT also requires that all models use the same data set:

```{r error=TRUE}
anova(fit_sex, fit_sex_toeph_subset, fit_sex_us)
```

#### Data sameness assertion: Different data, `refit = TRUE`

As the error message above suggests, we can set `refit = TRUE`, which will find the largest possible data set in common between all the fits, refit all the models using this data set, and then attempt to perform LRT:

```{r}
anova(fit_sex, fit_sex_toeph_subset, fit_sex_us, refit = TRUE)
```

We have success, and we see that the diagnostics have changed for Models 1 and 3 since they have been refitted with smaller data sets. Notice also the change in *df* for Model 3.

Concerning Model 1, the *df* are still 6 because the number of variance parameters in the `ar1` covariance structure is always 2. Therefore, changing the number of time points will not change the *df*; this also applies to `cs`, whose number of covariance parameters is also fixed at 2.

Concerning Model 2, the unstructured covariance has $m(m+1)/2$ parameters. In this case, when we account for the four coefficients, $df = 4 + 3(3+1)/2 = 10$.

Since Model 2's data set was a subset of the other Models' data, its results are unchanged.

#### Refitting without LRT

We can choose to refit the models while bypassing LRT results by specifying both `test = FALSE` and `refit = TRUE`. When we do this, there is no data sameness assertion.

```{r}
anova(fit_sex_toeph_subset, fit_sex_toeph, fit_race, fit_sex_csh,
      refit = TRUE,
      test = FALSE)
```

If we compare the last three models' results to the corresponding results above, we find different results because they were refitted with Model 1's smaller data set.

### Greater restrictions when using REML estimation

All model fits above utilize maximum likelihood (ML) estimation, and many of the above comparisons are only possible because no restricted maximum likelihood (REML) models were present. Specifically, if at least one model passed to `anova.mmrm()` used REML estimation, all models must have the same covariates. Covariance structures may still be nested.

#### Setup

Let's first create a few new models using we created above using REML by setting `reml = TRUE`:

```{r}
fit_sex_ar1_reml <- 
  mmrm(
    FEV1 ~ FEV1_BL + SEX + ARMCD + ar1(AVISIT | USUBJID),
    data = fev_data,
    reml = TRUE
  )
```

#### Bypassing LRT

REML models are not subject to the greater restrictions if LRT is bypassed.

As above, we can get results without refitting...

```{r}
anova(fit_sex_ar1_reml, fit_sex_toeph_subset, fit_interaction, test = FALSE)
```

...and with refitting:

```{r}
anova(fit_sex_ar1_reml, fit_sex_toeph_subset, fit_interaction,
      test = FALSE, refit = TRUE)
```

#### LRT with REML-estimated models

If we want LRT results with even one REML model, we will not be able to do so if models differ at all in their sets of covariates:

```{r error=TRUE}
anova(fit_sex_ar1_reml, fit_interaction)
```

But we can get LRT results if we are careful only to use models with the same covariates:

```{r}
anova(fit_sex_ar1_reml, fit_sex_us)
```

# Appendix: nested covariance structures

For all covariance structures described below, let $i$ and $j$ be the indices of the time points of a model. Therefore, $i$ and $j$ can each take a value in the sequence $\{1, 2, ...\}$, and the $i,j$ element of the covariance matrix is the variance between the $i$-th observation and the $j$-th observation.

## Unstructured covariance

Tautologically, all covariance structures are special cases of an unstructured covariance, and a model with a covariance structure can be considered "nested" within an model without a covariance structure, assuming that the covariates are also nested.

## Homogeneous and heterogeneous

All homogeneous covariance structures are special cases of their heterogeneous counterparts. Let $\sigma_i = \sigma$ and $\sigma_j = \sigma$

Therefore:

### `ad` is nested within `adh`:

$$
\sigma_{i}\sigma_{j}\prod_{k=i}^{j-1}\rho_{k} = \sigma\cdot\sigma\prod_{k=i}^{j-1}\rho_{k} = \sigma^{2}\prod_{k=i}^{j-1}\rho_{k}
$$

### `ar1` is nested within `ar1h`

$$
\sigma_{i}\sigma_{j}\rho^{\left \vert {i-j} \right \vert} = \sigma\cdot\sigma\rho^{\left \vert {i-j} \right \vert} = \sigma^{2}\rho^{\left \vert {i-j} \right \vert}
$$

### `cs` is nested within `csh`

$$
\sigma_{i}\sigma_{j}\left[ \rho I(i \neq j)+I(i=j) \right] = \sigma\cdot\sigma\left[ \rho I(i \neq j)+I(i=j) \right] = \sigma^{2}\left[ \rho I(i \neq j)+I(i=j) \right]
$$

### `toep` is nested within `toeph`

$$
\sigma_{i}\sigma_{j}\rho_{\left \vert {i-j} \right \vert} = \sigma\cdot\sigma\rho_{\left \vert {i-j} \right \vert} = \sigma^2\rho_{\left \vert {i-j} \right \vert}
$$

## Other hierarchies

We will now restrict the correlations in one covariance structure to produce another covariance structure in this set, demonstrating that they are nested. We will only do this for the homogeneous structures, since the heterogeneous ones plainly follow the same pattern.

### `ar1` is nested within `ad`

Let $\rho_k$ = $\rho$

Therefore:

$$
\sigma^2\prod_{k=i}^{j-1}\rho_{k} = \sigma^2\prod_{k=i}^{j-1}\rho = \sigma^2\rho^{\vert (j-1)-i+1\vert} = \sigma^2\rho^{\vert j-i\vert} = \sigma^2\rho^{\vert i-j\vert} 
$$

### `ar1` is nested within `toep`

Let $\rho$ be a value not necessarily in the sequence $\{\rho_0, \rho_1, \rho_2, ...\}$

Let $\rho_k = \rho^k$

Therefore:

$$
\sigma^2\rho_{\left \vert {i-j} \right \vert} = \sigma^{2}\rho^{\left \vert {i-j} \right \vert}
$$

### `cs` is nested within `toep`

Let $\rho$ be a value not necessarily in the sequence $\{\rho_0, \rho_1, \rho_2, ...\}$

Let $\rho_0 = 1$

Let $\rho_k=\rho$ $\forall$ $k \in \{1, 2, ...\}$

Observe that:

$$
x = x(I(true) + I(false)) = x\cdot I(true) + x\cdot I(false))
$$

Observe also that when $i \neq j$, $\left \vert i - j \right \vert \in \{1, 2, ...\}$, and when $i=j$, $\left \vert i - j \right \vert = 0$.

Therefore:

$$
\sigma^2\rho_{\left \vert {i-j} \right \vert} = \sigma^2\left[\rho_{\left \vert {i-j} \right \vert} I(i \neq j) + \rho_{\left \vert {i-j} \right \vert} I(i=j) \right] 
\\= \sigma^2\left[\rho I(i \neq j) + \rho_0 I(i=j) \right] 
\\= \sigma^2\left[\rho I(i \neq j) + 1 \cdot I(i=j) \right] 
\\= \sigma^{2}\left[ \rho I(i \neq j)+I(i=j) \right]
$$

## Resulting `list` object describing the covariance structure nesting hierarchy

```{r eval=FALSE}
list(
  ad = c("adh", "us"),
  adh = "us",
  ar1 = c("ad", "adh", "ar1h", "toep", "toeph", "us"),
  ar1h = c("adh", "toeph", "us"),
  cs = c("csh", "toep", "toeph", "us"),
  csh = c("toeph", "us"),
  toep = c("toeph", "us"),
  toeph = "us",
  us = character(0L),
  sp_exp = character(0L)
)
```

Note that since `sp_exp` is a valid covariance structure in `mmrm`, `anova.mmrm()` could encounter it. Therefore, it is included in this `list` object.
