mmrm coverage - 96.92%

Files
Source

#' Fitting an MMRM with Single Optimizer
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This function helps to fit an MMRM using `TMB` with a single optimizer,
#' while capturing messages and warnings.
#'
#' @inheritParams mmrm
#' @param control (`mmrm_control`)\cr object.
#' @param tmb_data (`mmrm_tmb_data`)\cr object.
#' @param formula_parts (`mmrm_tmb_formula_parts`)\cr object.
#' @param ... Additional arguments to pass to [mmrm_control()].
#'
#' @details
#' `fit_single_optimizer` will fit the `mmrm` model using the `control` provided.
#' If there are multiple optimizers provided in `control`, only the first optimizer
#' will be used.
#' If `tmb_data` and `formula_parts` are both provided, `formula`, `data`, `weights`,
#' `reml`, and `covariance` are ignored.
#'
#' @return The `mmrm_fit` object, with additional attributes containing warnings,
#'   messages, optimizer used and convergence status in addition to the
#'   `mmrm_tmb` contents.
#' @export
#'
#' @examples
#' mod_fit <- fit_single_optimizer(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data,
#'   weights = rep(1, nrow(fev_data)),
#'   optimizer = "nlminb"
#' )
#' attr(mod_fit, "converged")
fit_single_optimizer <- function(formula,
                                 data,
                                 weights,
                                 reml = TRUE,
                                 covariance = NULL,
                                 tmb_data,
                                 formula_parts,
                                 ...,
                                 control = mmrm_control(...)) {
  to_remove <- list(
    # Transient visit to invalid parameters.
    warnings = c("NA/NaN function evaluation")
  )
  as_diverged <- list(
    errors = c(
      "NA/NaN Hessian evaluation",
      "L-BFGS-B needs finite values of 'fn'"
    )
  )
  if (missing(tmb_data) || missing(formula_parts)) {
    h_valid_formula(formula)
    assert_data_frame(data)
    assert_numeric(weights, any.missing = FALSE, lower = .Machine$double.xmin)
    assert_flag(reml)
    assert_class(control, "mmrm_control")
    assert_list(control$optimizers, names = "unique", types = c("function", "partial"))
    quiet_fit <- h_record_all_output(
      fit_mmrm(
        formula = formula,
        data = data,
        weights = weights,
        reml = reml,
        covariance = covariance,
        control = control
      ),
      remove = to_remove,
      divergence = as_diverged
    )
  } else {
    assert_class(tmb_data, "mmrm_tmb_data")
    assert_class(formula_parts, "mmrm_tmb_formula_parts")
    quiet_fit <- h_record_all_output(
      fit_mmrm(
        formula_parts = formula_parts,
        tmb_data = tmb_data,
        control = control
      ),
      remove = to_remove,
      divergence = as_diverged
    )
  }
  if (length(quiet_fit$errors)) {
    stop(quiet_fit$errors)
  }
  converged <- (length(quiet_fit$warnings) == 0L) &&
    (length(quiet_fit$divergence) == 0L) &&
    isTRUE(quiet_fit$result$opt_details$convergence == 0)
  structure(
    quiet_fit$result,
    warnings = quiet_fit$warnings,
    messages = quiet_fit$messages,
    divergence = quiet_fit$divergence,
    converged = converged,
    class = c("mmrm_fit", class(quiet_fit$result))
  )
}

#' Summarizing List of Fits
#'
#' @param all_fits (`list` of `mmrm_fit` or `try-error`)\cr list of fits.
#'
#' @return List with `warnings`, `messages`, `log_liks` and `converged` results.
#' @keywords internal
h_summarize_all_fits <- function(all_fits) {
  assert_list(all_fits, types = c("mmrm_fit", "try-error"))
  is_error <- vapply(all_fits, is, logical(1), class2 = "try-error")

  warnings <- messages <- vector(mode = "list", length = length(all_fits))
  warnings[is_error] <- lapply(all_fits[is_error], as.character)
  warnings[!is_error] <- lapply(all_fits[!is_error], attr, which = "warnings")
  messages[!is_error] <- lapply(all_fits[!is_error], attr, which = "messages")
  log_liks <- as.numeric(rep(NA, length.out = length(all_fits)))
  log_liks[!is_error] <- vapply(all_fits[!is_error], stats::logLik, numeric(1L))
  converged <- rep(FALSE, length.out = length(all_fits))
  converged[!is_error] <- vapply(all_fits[!is_error], attr, logical(1), which = "converged")

  list(
    warnings = warnings,
    messages = messages,
    log_liks = log_liks,
    converged = converged
  )
}

#' Refitting MMRM with Multiple Optimizers
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @param fit (`mmrm_fit`)\cr original model fit from [fit_single_optimizer()].
#' @param ... Additional arguments passed to [mmrm_control()].
#' @param control (`mmrm_control`)\cr object.
#'
#' @return The best (in terms of log likelihood) fit which converged.
#'
#' @note For Windows, no parallel computations are currently implemented.
#' @export
#'
#' @examples
#' fit <- fit_single_optimizer(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data,
#'   weights = rep(1, nrow(fev_data)),
#'   optimizer = "nlminb"
#' )
#' best_fit <- refit_multiple_optimizers(fit)
refit_multiple_optimizers <- function(fit,
                                      ...,
                                      control = mmrm_control(...)) {
  assert_class(fit, "mmrm_fit")
  assert_class(control, "mmrm_control")

  n_cores_used <- ifelse(
    .Platform$OS.type == "windows",
    1L,
    min(
      length(control$optimizers),
      control$n_cores
    )
  )
  controls <- h_split_control(
    control,
    start = fit$theta_est
  )

  # Take the results from old fit as starting values for new fits.
  all_fits <- suppressWarnings(parallel::mcmapply(
    FUN = fit_single_optimizer,
    control = controls,
    MoreArgs = list(
      tmb_data = fit$tmb_data,
      formula_parts = fit$formula_parts
    ),
    mc.cores = n_cores_used,
    mc.silent = TRUE,
    SIMPLIFY = FALSE
  ))
  all_fits <- c(all_fits, list(old_result = fit))

  # Find the results that are ok and return best in terms of log-likelihood.
  all_fits_summary <- h_summarize_all_fits(all_fits)
  is_ok <- all_fits_summary$converged
  if (!any(is_ok)) {
    stop(
      "No optimizer led to a successful model fit. ",
      "Please try to use a different covariance structure or other covariates."
    )
  }
  best_optimizer <- which.max(all_fits_summary$log_liks[is_ok])
  all_fits[[which(is_ok)[best_optimizer]]]
}

#' Control Parameters for Fitting an MMRM
#'
#' @description `r lifecycle::badge("stable")`
#' Fine-grained specification of the MMRM fit details is possible using this
#' control function.
#'
#' @param n_cores (`count`)\cr number of cores to be used.
#' @param method (`string`)\cr adjustment method for degrees of freedom.
#' @param vcov (`string`)\cr coefficients covariance matrix adjustment method.
#' @param start (`NULL`, `numeric` or `function`)\cr optional start values for variance
#'   parameters. See details for more information.
#' @param accept_singular (`flag`)\cr whether singular design matrices are reduced
#'   to full rank automatically and additional coefficient estimates will be missing.
#' @param optimizers (`list`)\cr optimizer specification, created with [h_get_optimizers()].
#' @param drop_visit_levels (`flag`)\cr whether to drop levels for visit variable,
#'   if visit variable is a factor, see details.
#' @param ... additional arguments passed to [h_get_optimizers()].
#'
#' @details
#  - The `drop_visit_levels` flag will decide whether unobserved visits will be kept for analysis.
#'   For example, if the data only has observations at visits `VIS1`, `VIS3` and `VIS4`, by default
#'   they are treated to be equally spaced, the distance from `VIS1` to `VIS3`, and from `VIS3` to `VIS4`,
#'   are identical. However, you can manually convert this visit into a factor, with
#'   `levels = c("VIS1", "VIS2", "VIS3", "VIS4")`, and also use `drop_visits_levels = FALSE`,
#'   then the distance from `VIS1` to `VIS3` will be double, as `VIS2` is a valid visit.
#'   However, please be cautious because this can lead to convergence failure
#'   when using an unstructured covariance matrix and there are no observations
#'   at the missing visits.
#' - The `method` and `vcov` arguments specify the degrees of freedom and coefficients
#'   covariance matrix adjustment methods, respectively.
#'   - Allowed `vcov` includes: "Asymptotic", "Kenward-Roger", "Kenward-Roger-Linear", "Empirical" (CR0),
#'     "Empirical-Jackknife" (CR3), and "Empirical-Bias-Reduced" (CR2).
#'   - Allowed `method` includes: "Satterthwaite", "Kenward-Roger", "Between-Within" and "Residual".
#'   - If `method` is "Kenward-Roger" then only "Kenward-Roger" or "Kenward-Roger-Linear" are allowed for `vcov`.
#' - The `vcov` argument can be `NULL` to use the default covariance method depending on the `method`
#'   used for degrees of freedom, see the following table:
#'
#'    | `method`  |  Default `vcov`|
#'    |-----------|----------|
#'    |Satterthwaite| Asymptotic|
#'    |Kenward-Roger| Kenward-Roger|
#'    |Residual| Empirical|
#'    |Between-Within| Asymptotic|
#'
#' - Please note that "Kenward-Roger" for "Unstructured" covariance gives different results
#'   compared to SAS; Use "Kenward-Roger-Linear" for `vcov` instead for better matching
#'   of the SAS results.
#'
#' - The argument `start` is used to facilitate the choice of initial values for fitting the model.
#'   If `function` is provided, make sure its parameter is a valid element of `mmrm_tmb_data`
#'   or `mmrm_tmb_formula_parts` and it returns a numeric vector.
#'   By default or if `NULL` is provided, `std_start` will be used.
#'   Other implemented methods include `emp_start`.
#'
#' @return List of class `mmrm_control` with the control parameters.
#' @export
#'
#' @examples
#' mmrm_control(
#'   optimizer_fun = stats::optim,
#'   optimizer_args = list(method = "L-BFGS-B")
#' )
mmrm_control <- function(n_cores = 1L,
                         method = c("Satterthwaite", "Kenward-Roger", "Residual", "Between-Within"),
                         vcov = NULL,
                         start = std_start,
                         accept_singular = TRUE,
                         drop_visit_levels = TRUE,
                         ...,
                         optimizers = h_get_optimizers(...)) {
  assert_count(n_cores, positive = TRUE)
  assert_character(method)
  if (is.null(start)) {
    start <- std_start
  }
  assert(
    check_function(start, args = "..."),
    check_numeric(start, null.ok = FALSE),
    combine = "or"
  )
  assert_flag(accept_singular)
  assert_flag(drop_visit_levels)
  assert_list(optimizers, names = "unique", types = c("function", "partial"))
  assert_string(vcov, null.ok = TRUE)
  method <- match.arg(method)
  if (is.null(vcov)) {
    vcov <- h_get_cov_default(method)
  }
  assert_subset(
    vcov,
    c(
      "Asymptotic",
      "Empirical",
      "Empirical-Bias-Reduced",
      "Empirical-Jackknife",
      "Kenward-Roger",
      "Kenward-Roger-Linear"
    )
  )
  if (xor(identical(method, "Kenward-Roger"), vcov %in% c("Kenward-Roger", "Kenward-Roger-Linear"))) {
    stop(paste(
      "Kenward-Roger degrees of freedom must work together with Kenward-Roger",
      "or Kenward-Roger-Linear covariance!"
    ))
  }
  structure(
    list(
      optimizers = optimizers,
      start = start,
      accept_singular = accept_singular,
      method = method,
      vcov = vcov,
      n_cores = as.integer(n_cores),
      drop_visit_levels = drop_visit_levels
    ),
    class = "mmrm_control"
  )
}

#' Fit an MMRM
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This is the main function fitting the MMRM.
#'
#' @param formula (`formula`)\cr the model formula, see details.
#' @param data (`data`)\cr the data to be used for the model.
#' @param weights (`vector`)\cr an optional vector of weights to be used in
#'   the fitting process. Should be `NULL` or a numeric vector.
#' @param reml (`flag`)\cr whether restricted maximum likelihood (REML)
#'   estimation is used, otherwise maximum likelihood (ML) is used.
#' @param covariance (`cov_struct`)\cr a covariance structure type definition
#'   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 control (`mmrm_control`)\cr fine-grained fitting specifications list
#'   created with [mmrm_control()].
#' @param ... arguments passed to [mmrm_control()].
#'
#' @details
#' The `formula` typically looks like:
#' `FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID)`
#' so specifies response and covariates as usual, and exactly one special term
#' defines which covariance structure is used and what are the time point and
#' subject variables. The covariance structures in the formula can be
#' found in [`covariance_types`].
#'
#' The time points have to be unique for each subject. That is,
#' there cannot be time points with multiple observations for any subject.
#' The rationale is that these observations would need to be correlated, but it
#' is not possible within the currently implemented covariance structure framework
#' to do that correctly. Moreover, for non-spatial covariance structures, the time
#' variable must be a factor variable.
#'
#' When optimizer is not set, first the default optimizer
#' (`L-BFGS-B`) is used to fit the model. If that converges, this is returned.
#' If not, the other available optimizers from [h_get_optimizers()],
#' including `BFGS`, `CG` and `nlminb` are
#' tried (in parallel if `n_cores` is set and not on Windows).
#' If none of the optimizers converge, then the function fails. Otherwise
#' the best fit is returned.
#'
#' Note that fine-grained control specifications can either be passed directly
#' to the `mmrm` function, or via the `control` argument for bundling together
#' with the [mmrm_control()] function. Both cannot be used together, since
#' this would delete the arguments passed via `mmrm`.
#'
#' @return An `mmrm` object.
#'
#' @note The `mmrm` object is also an `mmrm_fit` and an `mmrm_tmb` object,
#' therefore corresponding methods also work (see [`mmrm_tmb_methods`]).
#'
#' Additional contents depend on the choice of the adjustment `method`:
#' - If Satterthwaite adjustment is used, the Jacobian information `jac_list`
#' is included.
#' - If Kenward-Roger adjustment is used, `kr_comp` contains necessary
#' components and `beta_vcov_adj` includes the adjusted coefficients covariance
#' matrix.
#'
#' Use of the package `emmeans` is supported, see [`emmeans_support`].
#'
#' NA values are always omitted regardless of `na.action` setting.
#'
#' When the number of visit levels is large, it usually requires large memory to create the
#' covariance matrix. By default, the maximum allowed visit levels is 100, and if there are more
#' visit levels, a confirmation is needed if run interactively.
#' You can use `options(mmrm.max_visits = <target>)` to increase the maximum allowed number of visit
#' levels. In non-interactive sessions the confirmation is not raised and will directly give you an error if
#' the number of visit levels exceeds the maximum.
#'
#' @export
#'
#' @examples
#' fit <- mmrm(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data
#' )
#'
#' # Direct specification of control details:
#' fit <- mmrm(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data,
#'   weights = fev_data$WEIGHTS,
#'   method = "Kenward-Roger"
#' )
#'
#' # Alternative specification via control argument (but you cannot mix the
#' # two approaches):
#' fit <- mmrm(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data,
#'   control = mmrm_control(method = "Kenward-Roger")
#' )
mmrm <- function(formula,
                 data,
                 weights = NULL,
                 covariance = NULL,
                 reml = TRUE,
                 control = mmrm_control(...),
                 ...) {
  assert_false(!missing(control) && !missing(...))
  assert_class(control, "mmrm_control")
  assert_list(control$optimizers, min.len = 1)

  if (control$method %in% c("Kenward-Roger", "Kenward-Roger-Linear") && !reml) {
    stop("Kenward-Roger only works for REML")
  }
  h_valid_formula(formula)
  covariance <- h_reconcile_cov_struct(formula, covariance)
  formula_parts <- h_mmrm_tmb_formula_parts(formula, covariance)
  h_tmb_warn_non_deterministic()

  if (!missing(data)) {
    attr(data, which = "dataname") <- toString(match.call()$data)
  } else {
    # na.action set to na.pass to allow data to be full; will be futher trimmed later
    data <- model.frame(formula_parts$full_formula, na.action = "na.pass")
  }

  if (is.null(weights)) {
    weights <- rep(1, nrow(data))
  } else {
    attr(weights, which = "dataname") <- deparse(match.call()$weights)
  }
  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
  )
  fit <- structure("", class = "try-error")
  names_all_optimizers <- names(control$optimizers)
  while (is(fit, "try-error") && length(control$optimizers) > 0) {
    fit <- fit_single_optimizer(
      tmb_data = tmb_data,
      formula_parts = formula_parts,
      control = control
    )
    if (is(fit, "try-error")) {
      warning(paste0(
        "Divergence with optimizer ", names(control$optimizers[1L]), " due to problems: ",
        toString(attr(fit, "divergence"))
      ))
    }
    control$optimizers <- control$optimizers[-1]
  }
  if (!attr(fit, "converged")) {
    more_optimizers <- length(control$optimizers) >= 1L
    if (more_optimizers) {
      fit <- refit_multiple_optimizers(
        fit = fit,
        control = control
      )
    } else {
      all_problems <- unlist(
        attributes(fit)[c("errors", "warnings")],
        use.names = FALSE
      )
      stop(paste0(
        "Chosen optimizers '", toString(names_all_optimizers), "' led to problems during model fit:\n",
        paste(paste0(seq_along(all_problems), ") ", all_problems), collapse = ";\n"), "\n",
        "Consider trying multiple or different optimizers."
      ))
    }
  }
  fit_msg <- attr(fit, "messages")
  if (!is.null(fit_msg)) {
    message(paste(fit_msg, collapse = "\n"))
  }
  fit$call <- match.call()
  fit$call$formula <- formula
  fit$method <- control$method
  fit$vcov <- control$vcov
  if (control$vcov %in% c("Kenward-Roger", "Kenward-Roger-Linear")) {
    fit$kr_comp <- h_get_kr_comp(fit$tmb_data, fit$theta_est)
    fit$beta_vcov_adj <- h_var_adj(
      v = fit$beta_vcov,
      w = component(fit, "theta_vcov"),
      p = fit$kr_comp$P,
      q = fit$kr_comp$Q,
      r = fit$kr_comp$R,
      linear = (control$vcov == "Kenward-Roger-Linear")
    )
  } else if (control$vcov %in% c("Empirical", "Empirical-Bias-Reduced", "Empirical-Jackknife")) {
    empirical_comp <- h_get_empirical(
      fit$tmb_data, fit$theta_est, fit$beta_est, fit$beta_vcov, control$vcov
    )
    fit$beta_vcov_adj <- empirical_comp$cov
    fit$empirical_df_mat <- empirical_comp$df_mat
    fit$score_per_subject <- empirical_comp$score_per_subject
    dimnames(fit$beta_vcov_adj) <- dimnames(fit$beta_vcov)
  } else if (identical(control$vcov, "Asymptotic")) {
    # Note that we only need the Jacobian list under Asymptotic covariance method,
    # cf. the Satterthwaite vignette.
    if (identical(fit$method, "Satterthwaite")) {
      fit$jac_list <- h_jac_list(fit$tmb_data, fit$theta_est, fit$beta_vcov)
    }
  } else {
    stop("Unrecognized coefficent variance-covariance method!")
  }

  class(fit) <- c("mmrm", class(fit))
  fit
}

#' Register `mmrm` For Use With `car::Anova`
#'
#' @inheritParams base::requireNamespace
#' @return A logical value indicating whether registration was successful.
#'
#' @keywords internal
car_add_mmrm <- function(quietly = FALSE) {
  if (!requireNamespace("car", quietly = quietly)) {
    return(FALSE)
  }
  envir <- asNamespace("mmrm")
  h_register_s3("car", "Anova", "mmrm", envir)
  TRUE
}


#' Obtain Contrast for Specified Effect
#'
#' This is support function to obtain contrast matrix for type II/III testing.
#'
#' @param object (`mmrm`)\cr the fitted MMRM.
#' @param effect (`string`) the name of the effect.
#' @param type (`string`) type of test, "II", "III", '2', or '3'.
#' @param tol (`numeric`) threshold blow which values are treated as 0.
#'
#' @return A `matrix` of the contrast.
#'
#' @keywords internal
h_get_contrast <- function(object, effect, type = c("II", "III", "2", "3"), tol = sqrt(.Machine$double.eps)) {
  assert_class(object, "mmrm")
  assert_string(effect)
  assert_double(tol, finite = TRUE, len = 1L)
  type <- match.arg(type)
  mx <- component(object, "x_matrix")
  asg <- attr(mx, "assign")
  formula <- object$formula_parts$model_formula
  tms <- terms(formula)
  fcts <- attr(tms, "factors")[-1L, , drop = FALSE] # Discard the response.
  ods <- attr(tms, "order")
  assert_subset(effect, colnames(fcts))
  idx <- which(effect == colnames(fcts))
  cols <- which(asg == idx)
  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)
  }
  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)) {
    additional_vars <- names(which(fcts[, i] > fcts[, idx]))
    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)
          ret <- solve(t(x1) %*% m %*% x1) %*% t(x1) %*% m %*% x2
          if (contains_intercept) {
            ret[-1, ] - ret[1, ]
          } else {
            ret
          }
        },
        "3" = ,
        "III" = {
          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
    }
  }
  l_mx[abs(l_mx) < tol] <- 0
  l_mx
}

#' Conduct type II/III hypothesis testing on the MMRM fit results.
#'
#' @param mod (`mmrm`)\cr the fitted MMRM.
#' @param ... not used.
#' @inheritParams h_get_contrast
#'
#' @details
#' `Anova` will return `anova` object with one row per variable and columns
#' `Num Df`(numerator degrees of freedom), `Denom Df`(denominator degrees of freedom),
#' `F Statistic` and `Pr(>=F)`.
#'
#' @keywords internal
# Please do not load `car` and then create the documentation. The Rd file will be different.
Anova.mmrm <- function(mod, type = c("II", "III", "2", "3"), tol = sqrt(.Machine$double.eps), ...) { # nolint
  assert_double(tol, finite = TRUE, len = 1L)
  type <- match.arg(type)
  vars <- colnames(attr(terms(mod$formula_parts$model_formula), "factors"))
  ret <- lapply(
    vars,
    function(x) df_md(mod, h_get_contrast(mod, x, type, tol))
  )
  ret_df <- do.call(rbind.data.frame, ret)
  row.names(ret_df) <- vars
  colnames(ret_df) <- c("Num Df", "Denom Df", "F Statistic", "Pr(>=F)")
  class(ret_df) <- c("anova", "data.frame")
  attr(ret_df, "heading") <- sprintf(
    "Analysis of Fixed Effect Table (Type %s F tests)",
    switch(type,
      "2" = ,
      "II" = "II",
      "3" = ,
      "III" = "III"
    )
  )
  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
  )
}

#' Tidying Methods for `mmrm` Objects
#'
#' @description `r lifecycle::badge("stable")`
#'
#' These methods tidy the estimates from an `mmrm` object into a
#' summary.
#'
#' @param x (`mmrm`)\cr fitted model.
#' @param conf.int (`flag`)\cr if `TRUE` columns for the lower (`conf.low`) and upper bounds
#'   (`conf.high`) of coefficient estimates are included.
#' @param conf.level (`number`)\cr defines the range of the optional confidence internal.
#' @param newdata (`data.frame` or `NULL`)\cr optional new data frame.
#' @param se_fit (`flag`)\cr whether to return standard errors of fit.
#' @param interval (`string`)\cr type of interval calculation.
#' @param type.residuals (`string`)\cr passed on to [residuals.mmrm_tmb()].
#' @param ... only used by `augment()` to pass arguments to the [predict.mmrm_tmb()] method.
#'
#' @name mmrm_tidiers
#' @aliases mmrm_tidiers
#'
#' @seealso [`mmrm_methods`], [`mmrm_tmb_methods`] for additional methods.
#'
#' @examples
#' fit <- mmrm(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data
#' )
NULL

#' @describeIn mmrm_tidiers derives tidy `tibble` from an `mmrm` object.
#' @exportS3Method
#' @examples
#' # Applying tidy method to return summary table of covariate estimates.
#' fit |> tidy()
#' fit |> tidy(conf.int = TRUE, conf.level = 0.9)
tidy.mmrm <- function(x, # nolint
                      conf.int = FALSE, # nolint
                      conf.level = 0.95, # nolint
                      ...) {
  assert_flag(conf.int)
  assert_number(conf.level, lower = 0, upper = 1)
  tbl <- tibble::as_tibble(summary(x)$coefficients, rownames = "term")
  colnames(tbl) <- c("term", "estimate", "std.error", "df", "statistic", "p.value")
  coefs <- coef(x)
  if (length(coefs) != nrow(tbl)) {
    coefs <- tibble::enframe(coefs, name = "term", value = "estimate")
    tbl <- merge(coefs, tbl, by = c("term", "estimate"))
  }
  if (conf.int) {
    ci <- h_tbl_confint_terms(x, level = conf.level)
    tbl <- tibble::as_tibble(merge(tbl, ci, by = "term"))
  }
  tbl
}

#' @describeIn mmrm_tidiers derives `glance` `tibble` from an `mmrm` object.
#' @exportS3Method
#' @examples
#' # Applying glance method to return summary table of goodness of fit statistics.
#' fit |> glance()
glance.mmrm <- function(x, ...) { # nolint
  tibble::as_tibble(summary(x)$aic_list)
}

#' @describeIn mmrm_tidiers derives `augment` `tibble` from an `mmrm` object.
#' @exportS3Method
#' @examples
#' # Applying augment method to return merged `tibble` of model data, fitted and residuals.
#' fit |> augment()
#' fit |> augment(interval = "confidence")
#' fit |> augment(type.residuals = "pearson")
augment.mmrm <- function(x, # nolint
                         newdata = NULL,
                         interval = c("none", "confidence", "prediction"),
                         se_fit = (interval != "none"),
                         type.residuals = c("response", "pearson", "normalized"), # nolint
                         ...) {
  type.residuals <- match.arg(type.residuals) # nolint
  resid_df <- NULL
  if (is.null(newdata)) {
    newdata <- stats::get_all_vars(x, data = stats::na.omit(x$data))
    resid_df <- data.frame(
      .rownames = rownames(newdata),
      .resid = unname(residuals(x, type = type.residuals))
    )
  }
  interval <- match.arg(interval)

  tbl <- h_newdata_add_pred(
    x,
    newdata = newdata,
    se_fit = se_fit,
    interval = interval,
    ...
  )
  if (!is.null(resid_df)) {
    tbl <- merge(tbl, resid_df, by = ".rownames")
    tbl$.rownames <- as.numeric(tbl$.rownames)
    tbl <- tbl[order(tbl$.rownames), , drop = FALSE]
  }
  tibble::as_tibble(tbl)
}

#' Extract `tibble` with Confidence Intervals and Term Names
#'
#' This is used in [tidy.mmrm()].
#'
#' @param x (`mmrm`)\cr fit object.
#' @param ... passed to [stats::confint()], hence not used at the moment.
#'
#' @return A `tibble` with `term`, `conf.low`, `conf.high` columns.
#'
#' @keywords internal
h_tbl_confint_terms <- function(x, ...) {
  df <- stats::confint(x, ...)
  tbl <- tibble::as_tibble(df, rownames = "term", .name_repair = "minimal")
  names(tbl) <- c("term", "conf.low", "conf.high")
  tbl
}

#' Add Prediction Results to New Data
#'
#' This is used in [augment.mmrm()].
#'
#' @param x (`mmrm`)\cr fit.
#' @param newdata (`data.frame`)\cr data to predict.
#' @param se_fit (`flag`)\cr whether to return standard error of prediction,
#'   can only be used when `interval` is not "none".
#' @param interval (`string`)\cr type of interval.
#' @param ... passed to [predict.mmrm_tmb()].
#'
#' @return The `newdata` as a `tibble` with additional columns `.fitted`,
#'   `.lower`, `.upper` (if interval is not `none`) and `.se.fit` (if `se_fit`
#'   requested).
#'
#' @keywords internal
h_newdata_add_pred <- function(x,
                               newdata,
                               se_fit,
                               interval,
                               ...) {
  assert_class(x, "mmrm")
  assert_data_frame(newdata)
  assert_flag(se_fit)
  assert_string(interval)
  if (interval == "none") {
    assert_false(se_fit)
  }

  tbl <- h_df_to_tibble(newdata)
  pred_results <- predict(
    x,
    newdata = newdata,
    na.action = stats::na.pass,
    se.fit = se_fit,
    interval = interval,
    ...
  )
  if (interval == "none") {
    assert_numeric(pred_results)
    tbl$.fitted <- unname(pred_results)
  } else {
    assert_matrix(pred_results)
    tbl$.fitted <- unname(pred_results[, "fit"])
    tbl$.lower <- unname(pred_results[, "lwr"])
    tbl$.upper <- unname(pred_results[, "upr"])
  }
  if (se_fit) {
    tbl$.se.fit <- unname(pred_results[, "se"])
  }
  tbl
}

#' Coerce a Data Frame to a `tibble`
#'
#' This is used in [h_newdata_add_pred()].
#'
#' @details This is only a thin wrapper around [tibble::as_tibble()], except
#' giving a useful error message and it checks for `rownames` and adds them
#' as a new column `.rownames` if they are not just a numeric sequence as
#' per the [tibble::has_rownames()] decision.
#'
#' @param data (`data.frame`)\cr what to coerce.
#'
#' @return The `data` as a `tibble`, potentially with a `.rownames` column.
#'
#' @keywords internal
h_df_to_tibble <- function(data) {
  tryCatch(tbl <- tibble::as_tibble(data), error = function(cnd) {
    stop("Could not coerce data to `tibble`. Try explicitly passing a",
      "dataset to either the `data` or `newdata` argument.",
      call. = FALSE
    )
  })
  if (tibble::has_rownames(data)) {
    tbl <- tibble::add_column(tbl, .rownames = rownames(data), .before = TRUE)
  }
  tbl
}

#' Processing the Formula for `TMB` Fit
#'
#' @param formula (`formula`)\cr Original formula.
#' @param covariance (`cov_struct`)\cr A covariance structure from which
#'   additional formula parts should be added.
#'
#' @return List of class `mmrm_tmb_formula_parts` with elements:
#'
#' - `formula`: the original input.
#' - `model_formula`: `formula` with the covariance term is removed.
#' - `model_formula`: `formula` with the covariance term removed.
#' - `full_formula`: same as `model_formula` but includes the covariance
#'   structure's subject, visit and (optionally) group variables.
#' - `cov_type`: `string` with covariance term type (e.g. `"us"`).
#' - `is_spatial`: `flag` indicator of whether the covariance structure is
#'   spatial
#' - `visit_var`: `character` with the visit variable name.
#' - `subject_var`: `string` with the subject variable name.
#' - `group_var`: `string` with the group variable name. If no group specified,
#'   this element is `NULL`.
#' - `model_var`: `character` with the variables names of the formula, except `subject_var`.
#'
#' @keywords internal
h_mmrm_tmb_formula_parts <- function(
    formula,
    covariance = as.cov_struct(formula, warn_partial = FALSE)) {
  assert_formula(formula)
  assert_true(identical(length(formula), 3L))

  model_formula <- h_drop_covariance_terms(formula)

  structure(
    list(
      formula = formula,
      model_formula = model_formula,
      full_formula = h_add_covariance_terms(model_formula, covariance),
      cov_type = tmb_cov_type(covariance),
      is_spatial = covariance$type == "sp_exp",
      visit_var = covariance$visits,
      subject_var = covariance$subject,
      group_var = if (length(covariance$group) < 1) NULL else covariance$group,
      model_var = setdiff(all.vars(formula[[3]]), covariance$subject)
    ),
    class = "mmrm_tmb_formula_parts"
  )
}

#' Data for `TMB` Fit
#'
#' @param formula_parts (`mmrm_tmb_formula_parts`)\cr list with formula parts
#'   from [h_mmrm_tmb_formula_parts()].
#' @param data (`data.frame`)\cr which contains variables used in `formula_parts`.
#' @param weights (`vector`)\cr weights to be used in the fitting process.
#' @param reml (`flag`)\cr whether restricted maximum likelihood (REML) estimation is used,
#'   otherwise maximum likelihood (ML) is used.
#' @param singular (`string`)\cr choices of method deal with rank-deficient matrices. "error" to
#'   stop the function return the error, "drop" to drop these columns, and "keep" to keep all the columns.
#' @param drop_visit_levels (`flag`)\cr whether to drop levels for visit variable, if visit variable is a factor.
#' @param allow_na_response (`flag`)\cr whether NA in response is allowed.
#' @param drop_levels (`flag`)\cr whether drop levels for covariates. If not dropped could lead to singular matrix.
#'
#' @return List of class `mmrm_tmb_data` with elements:
#' - `full_frame`: `data.frame` with `n` rows containing all variables needed in the model.
#' - `data`: `data.frame` of input dataset.
#' - `x_matrix`: `matrix` with `n` rows and `p` columns specifying the overall design matrix.
#' - `x_cols_aliased`: `logical` with potentially more than `p` elements indicating which
#'      columns in the original design matrix have been left out to obtain a full rank
#'      `x_matrix`.
#' - `y_vector`: length `n` `numeric` specifying the overall response vector.
#' - `weights_vector`: length `n` `numeric` specifying the weights vector.
#' - `n_visits`: `int` with the number of visits, which is the dimension of the
#'      covariance matrix.
#' - `n_subjects`: `int` with the number of subjects.
#' - `subject_zero_inds`: length `n_subjects` `integer` containing the zero-based start
#'     indices for each subject.
#' - `subject_n_visits`: length `n_subjects` `integer` containing the number of
#'     observed visits for each subjects. So the sum of this vector equals `n`.
#' - `cov_type`: `string` value specifying the covariance type.
#' - `is_spatial_int`: `int` specifying whether the covariance structure is spatial(1) or not(0).
#' - `reml`: `int` specifying whether REML estimation is used (1), otherwise ML (0).
#' - `subject_groups`: `factor` specifying the grouping for each subject.
#' - `n_groups`: `int` with the number of total groups
#'
#' @details Note that the `subject_var` must not be factor but can also be character.
#'   If it is character, then it will be converted to factor internally. Here
#'   the levels will be the unique values, sorted alphabetically and numerically if there
#'   is a common string prefix of numbers in the character elements. For full control
#'   on the order please use a factor.
#'
#' @keywords internal
h_mmrm_tmb_data <- function(formula_parts,
                            data,
                            weights,
                            reml,
                            singular = c("drop", "error", "keep"),
                            drop_visit_levels,
                            allow_na_response = FALSE,
                            drop_levels = TRUE,
                            xlev = NULL,
                            contrasts = NULL) {
  assert_class(formula_parts, "mmrm_tmb_formula_parts")
  assert_data_frame(data)
  varname <- formula_parts[grepl("_var", names(formula_parts))]
  assert_names(
    names(data),
    must.include = unlist(varname, use.names = FALSE)
  )
  assert_true(is.factor(data[[formula_parts$subject_var]]) || is.character(data[[formula_parts$subject_var]]))
  assert_numeric(weights, len = nrow(data))
  assert_flag(reml)
  singular <- match.arg(singular)
  assert_flag(drop_visit_levels)

  if (is.character(data[[formula_parts$subject_var]])) {
    data[[formula_parts$subject_var]] <- factor(
      data[[formula_parts$subject_var]],
      levels = stringr::str_sort(unique(data[[formula_parts$subject_var]]), numeric = TRUE)
    )
  }
  data_order <- if (formula_parts$is_spatial) {
    order(data[[formula_parts$subject_var]])
  } else {
    subject_visit_data <- data[, c(formula_parts$subject_var, formula_parts$visit_var)]
    is_duplicated <- duplicated(subject_visit_data)
    if (any(is_duplicated)) {
      stop(
        "time points have to be unique for each subject, detected following duplicates in data:\n",
        paste(utils::capture.output(print(subject_visit_data[is_duplicated, ])), collapse = "\n")
      )
    }
    order(data[[formula_parts$subject_var]], data[[formula_parts$visit_var]])
  }
  if (identical(formula_parts$is_spatial, FALSE)) {
    h_confirm_large_levels(length(levels(data[[formula_parts$visit_var]])))
  }
  data <- data[data_order, ]
  weights <- weights[data_order]
  data <- data.frame(data, weights)
  # Weights is always the last column.
  weights_name <- colnames(data)[ncol(data)]
  # If `y` is allowed to be NA, then first replace y with 1:n, then replace it with original y.
  if (!allow_na_response) {
    h_warn_na_action()
  }
  full_frame <- eval(
    bquote(stats::model.frame(
      formula_parts$full_formula,
      data = data,
      weights = .(as.symbol(weights_name)),
      na.action = "na.pass",
      xlev = xlev
    ))
  )
  if (drop_levels) {
    full_frame <- h_drop_levels(full_frame, formula_parts$subject_var, formula_parts$visit_var, names(xlev))
  }
  has_response <- !identical(attr(attr(full_frame, "terms"), "response"), 0L)
  keep_ind <- if (allow_na_response && has_response) {
    # Note that response is always the first column if there is response.
    stats::complete.cases(full_frame[, -1L, drop = FALSE])
  } else {
    stats::complete.cases(full_frame)
  }
  full_frame <- full_frame[keep_ind, ]
  if (drop_visit_levels && !formula_parts$is_spatial && h_extra_levels(full_frame[[formula_parts$visit_var]])) {
    visit_vec <- full_frame[[formula_parts$visit_var]]
    old_levels <- levels(visit_vec)
    full_frame[[formula_parts$visit_var]] <- droplevels(visit_vec)
    new_levels <- levels(full_frame[[formula_parts$visit_var]])
    dropped <- setdiff(old_levels, new_levels)
    message(
      "In ", formula_parts$visit_var, " there are dropped visits: ", toString(dropped),
      ".\n Additional attributes including contrasts are lost.\n",
      "To avoid this behavior, make sure use `drop_visit_levels = FALSE`."
    )
  }
  is_factor_col <- vapply(full_frame, is.factor, FUN.VALUE = TRUE)
  is_factor_col <- intersect(names(is_factor_col)[is_factor_col], all.vars(formula_parts$model_formula))
  x_matrix <- stats::model.matrix(
    formula_parts$model_formula,
    data = full_frame,
    contrasts.arg = h_default_value(contrasts, lapply(full_frame[is_factor_col], contrasts))
  )
  x_cols_aliased <- stats::setNames(rep(FALSE, ncol(x_matrix)), nm = colnames(x_matrix))
  qr_x_mat <- qr(x_matrix)
  if (qr_x_mat$rank < ncol(x_matrix)) {
    cols_to_drop <- utils::tail(qr_x_mat$pivot, ncol(x_matrix) - qr_x_mat$rank)
    if (identical(singular, "error")) {
      stop(
        "design matrix only has rank ", qr_x_mat$rank, " and ", length(cols_to_drop),
        " columns (", toString(colnames(x_matrix)[cols_to_drop]), ") could be dropped",
        " to achieve full rank ", ncol(x_matrix), " by using `accept_singular = TRUE`"
      )
    } else if (identical(singular, "drop")) {
      assign_attr <- attr(x_matrix, "assign")
      contrasts_attr <- attr(x_matrix, "contrasts")
      x_matrix <- x_matrix[, -cols_to_drop, drop = FALSE]
      x_cols_aliased[cols_to_drop] <- TRUE
      attr(x_matrix, "assign") <- assign_attr[-cols_to_drop]
      attr(x_matrix, "contrasts") <- contrasts_attr
    }
  }
  y_vector <- if (has_response) {
    as.numeric(stats::model.response(full_frame))
  } else {
    rep(NA_real_, nrow(full_frame))
  }
  weights_vector <- as.numeric(stats::model.weights(full_frame))
  n_subjects <- length(unique(full_frame[[formula_parts$subject_var]]))
  subject_zero_inds <- which(!duplicated(full_frame[[formula_parts$subject_var]])) - 1L
  subject_n_visits <- c(utils::tail(subject_zero_inds, -1L), nrow(full_frame)) - subject_zero_inds
  # It is possible that `subject_var` is factor with more levels (and this does not affect fit)
  # so no check is needed for `subject_visits`.
  assert_true(all(subject_n_visits > 0))
  if (!is.null(formula_parts$group_var)) {
    assert_factor(data[[formula_parts$group_var]])
    subject_groups <- full_frame[[formula_parts$group_var]][subject_zero_inds + 1L]
    n_groups <- nlevels(subject_groups)
  } else {
    subject_groups <- factor(rep(0L, n_subjects))
    n_groups <- 1L
  }
  coordinates <- full_frame[, formula_parts$visit_var, drop = FALSE]
  if (formula_parts$is_spatial) {
    lapply(coordinates, assert_numeric)
    coordinates_matrix <- as.matrix(coordinates)
    n_visits <- max(subject_n_visits)
  } else {
    assert(identical(ncol(coordinates), 1L))
    assert_factor(coordinates[[1L]])
    coordinates_matrix <- as.matrix(as.integer(coordinates[[1L]]) - 1, ncol = 1)
    n_visits <- nlevels(coordinates[[1L]])
    assert_true(all(subject_n_visits <= n_visits))
  }
  structure(
    list(
      full_frame = full_frame,
      data = data,
      x_matrix = x_matrix,
      x_cols_aliased = x_cols_aliased,
      coordinates = coordinates_matrix,
      y_vector = y_vector,
      weights_vector = weights_vector,
      n_visits = n_visits,
      n_subjects = n_subjects,
      subject_zero_inds = subject_zero_inds,
      subject_n_visits = subject_n_visits,
      cov_type = formula_parts$cov_type,
      is_spatial_int = as.integer(formula_parts$is_spatial),
      reml = as.integer(reml),
      subject_groups = subject_groups,
      n_groups = n_groups
    ),
    class = "mmrm_tmb_data"
  )
}

#' Start Parameters for `TMB` Fit
#'
#' @param formula_parts (`mmrm_tmb_formula_parts`)\cr produced by
#'  [h_mmrm_tmb_formula_parts()].
#' @param tmb_data (`mmrm_tmb_data`)\cr produced by [h_mmrm_tmb_data()].
#' @param start (`numeric` or `NULL`)\cr optional start values for variance
#'   parameters.
#' @param n_groups (`int`)\cr number of groups.
#' @return List with element `theta` containing the start values for the variance
#'   parameters.
#'
#' @keywords internal
h_mmrm_tmb_parameters <- function(formula_parts,
                                  tmb_data,
                                  start,
                                  n_groups = 1L) {
  assert_class(formula_parts, "mmrm_tmb_formula_parts")
  assert_class(tmb_data, "mmrm_tmb_data")

  m <- tmb_data$n_visits
  start_value0 <- std_start(formula_parts$cov_type, m, n_groups)
  theta_dim <- length(start_value0)
  start_values <- if (is.null(start)) {
    start_value0
  } else if (test_function(start)) {
    do.call(start, utils::modifyList(formula_parts, tmb_data))
  } else {
    start
  }
  assert_numeric(start_values, len = theta_dim, any.missing = FALSE, finite = TRUE)
  list(theta = start_values)
}

#' Asserting Sane Start Values for `TMB` Fit
#'
#' @param tmb_object (`list`)\cr created with [TMB::MakeADFun()].
#'
#' @return Nothing, only used for assertions.
#'
#' @keywords internal
h_mmrm_tmb_assert_start <- function(tmb_object) {
  assert_list(tmb_object)
  assert_subset(c("fn", "gr", "par"), names(tmb_object))

  if (is.na(tmb_object$fn(tmb_object$par))) {
    stop("negative log-likelihood is NaN at starting parameter values")
  }
  if (any(is.na(tmb_object$gr(tmb_object$par)))) {
    stop("some elements of gradient are NaN at starting parameter values")
  }
}

#' Checking the `TMB` Optimization Result
#'
#' @param tmb_opt (`list`)\cr optimization result.
#' @param mmrm_tmb (`mmrm_tmb`)\cr result from [h_mmrm_tmb_fit()].
#'
#' @return Nothing, only used to generate warnings in case that the model
#' did not converge.
#'
#' @keywords internal
h_mmrm_tmb_check_conv <- function(tmb_opt, mmrm_tmb) {
  assert_list(tmb_opt)
  assert_subset(c("par", "objective", "convergence", "message"), names(tmb_opt))
  assert_class(mmrm_tmb, "mmrm_tmb")

  if (!is.null(tmb_opt$convergence) && tmb_opt$convergence != 0) {
    warning("Model convergence problem: ", tmb_opt$message, ".")
    return()
  }
  theta_vcov <- mmrm_tmb$theta_vcov
  if (is(theta_vcov, "try-error")) {
    warning("Model convergence problem: hessian is singular, theta_vcov not available.")
    return()
  }
  if (!all(is.finite(theta_vcov))) {
    warning("Model convergence problem: theta_vcov contains non-finite values.")
    return()
  }
  eigen_vals <- eigen(theta_vcov, only.values = TRUE)$values
  if (mode(eigen_vals) == "complex" || any(eigen_vals <= 0)) {
    # Note: complex eigen values signal that the matrix is not symmetric, therefore not positive definite.
    warning("Model convergence problem: theta_vcov is not positive definite.")
    return()
  }
  qr_rank <- qr(theta_vcov)$rank
  if (qr_rank < ncol(theta_vcov)) {
    warning("Model convergence problem: theta_vcov is numerically singular.")
  }
}

#' Extract covariance matrix from `TMB` report and input data
#'
#' This helper does some simple post-processing to extract covariance matrix or named
#' list of covariance matrices if the fitting is using grouped covariance matrices.
#'
#' @param tmb_report (`list`)\cr report created with [TMB::MakeADFun()] report function.
#' @param tmb_data (`mmrm_tmb_data`)\cr produced by [h_mmrm_tmb_data()].
#' @param visit_var (`character`)\cr character vector of the visit variable
#' @param is_spatial (`flag`)\cr indicator whether the covariance structure is spatial.
#' @return Return a simple covariance matrix if there is no grouping, or a named
#' list of estimated grouped covariance matrices,
#' with its name equal to the group levels.
#'
#' @keywords internal
h_mmrm_tmb_extract_cov <- function(tmb_report, tmb_data, visit_var, is_spatial) {
  d <- dim(tmb_report$covariance_lower_chol)
  visit_names <- if (!is_spatial) {
    levels(tmb_data$full_frame[[visit_var]])
  } else {
    c(0, 1)
  }
  cov <- lapply(
    seq_len(d[1] / d[2]),
    function(i) {
      ret <- tcrossprod(tmb_report$covariance_lower_chol[seq(1 + (i - 1) * d[2], i * d[2]), ])
      dimnames(ret) <- list(visit_names, visit_names)
      return(ret)
    }
  )
  if (identical(tmb_data$n_groups, 1L)) {
    cov <- cov[[1]]
  } else {
    names(cov) <- levels(tmb_data$subject_groups)
  }
  return(cov)
}

#' Build `TMB` Fit Result List
#'
#' This helper does some simple post-processing of the `TMB` object and
#' optimization results, including setting names, inverting matrices etc.
#'
#' @param tmb_object (`list`)\cr created with [TMB::MakeADFun()].
#' @param tmb_opt (`list`)\cr optimization result.
#' @param formula_parts (`mmrm_tmb_formula_parts`)\cr produced by
#'  [h_mmrm_tmb_formula_parts()].
#' @param tmb_data (`mmrm_tmb_data`)\cr produced by [h_mmrm_tmb_data()].
#'
#' @return List of class `mmrm_tmb` with:
#'   - `cov`: estimated covariance matrix, or named list of estimated group specific covariance matrices.
#'   - `beta_est`: vector of coefficient estimates.
#'   - `beta_vcov`: Variance-covariance matrix for coefficient estimates.
#'   - `beta_vcov_inv_L`: Lower triangular matrix `L` of the inverse variance-covariance matrix decomposition.
#'   - `beta_vcov_inv_D`: vector of diagonal matrix `D` of the inverse variance-covariance matrix decomposition.
#'   - `theta_est`: vector of variance parameter estimates.
#'   - `theta_vcov`: variance-covariance matrix for variance parameter estimates.
#'   - `neg_log_lik`: obtained negative log-likelihood.
#'   - `formula_parts`: input.
#'   - `data`: input.
#'   - `weights`: input.
#'   - `reml`: input as a flag.
#'   - `opt_details`: list with optimization details including convergence code.
#'   - `tmb_object`: original `TMB` object created with [TMB::MakeADFun()].
#'   - `tmb_data`: input.
#'
#' @details Instead of inverting or decomposing `beta_vcov`, it can be more efficient to use its robust
#'   Cholesky decomposition `LDL^T`, therefore we return the corresponding two components `L` and `D`
#'   as well since they have been available on the `C++` side already.
#'
#' @keywords internal
h_mmrm_tmb_fit <- function(tmb_object,
                           tmb_opt,
                           formula_parts,
                           tmb_data) {
  assert_list(tmb_object)
  assert_subset(c("fn", "gr", "par", "he"), names(tmb_object))
  assert_list(tmb_opt)
  assert_subset(c("par", "objective", "convergence", "message"), names(tmb_opt))
  assert_class(formula_parts, "mmrm_tmb_formula_parts")
  assert_class(tmb_data, "mmrm_tmb_data")

  tmb_report <- tmb_object$report(par = tmb_opt$par)
  x_matrix_cols <- colnames(tmb_data$x_matrix)
  cov <- h_mmrm_tmb_extract_cov(tmb_report, tmb_data, formula_parts$visit_var, formula_parts$is_spatial)
  beta_est <- tmb_report$beta
  names(beta_est) <- x_matrix_cols
  beta_vcov <- tmb_report$beta_vcov
  dimnames(beta_vcov) <- list(x_matrix_cols, x_matrix_cols)
  beta_vcov_inv_L <- tmb_report$XtWX_L # nolint
  beta_vcov_inv_D <- tmb_report$XtWX_D # nolint
  theta_est <- tmb_opt$par
  names(theta_est) <- NULL
  theta_vcov <- try(solve(tmb_object$he(tmb_opt$par)), silent = TRUE)
  opt_details_names <- setdiff(
    names(tmb_opt),
    c("par", "objective")
  )
  structure(
    list(
      cov = cov,
      beta_est = beta_est,
      beta_vcov = beta_vcov,
      beta_vcov_inv_L = beta_vcov_inv_L,
      beta_vcov_inv_D = beta_vcov_inv_D,
      theta_est = theta_est,
      theta_vcov = theta_vcov,
      neg_log_lik = tmb_opt$objective,
      formula_parts = formula_parts,
      data = tmb_data$data,
      weights = tmb_data$weights_vector,
      reml = as.logical(tmb_data$reml),
      opt_details = tmb_opt[opt_details_names],
      tmb_object = tmb_object,
      tmb_data = tmb_data
    ),
    class = "mmrm_tmb"
  )
}

#' Low-Level Fitting Function for MMRM
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This is the low-level function to fit an MMRM. Note that this does not
#' try different optimizers or adds Jacobian information etc. in contrast to
#' [mmrm()].
#'
#' @param formula (`formula`)\cr model formula with exactly one special term
#'   specifying the visits within subjects, see details.
#' @param data (`data.frame`)\cr input data containing the variables used in
#'   `formula`.
#' @param weights (`vector`)\cr input vector containing the weights.
#' @inheritParams h_mmrm_tmb_data
#' @param covariance (`cov_struct`)\cr A covariance structure type definition,
#'   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 control (`mmrm_control`)\cr list of control options produced by
#'   [mmrm_control()].
#' @inheritParams fit_single_optimizer
#'
#' @return List of class `mmrm_tmb`, see [h_mmrm_tmb_fit()] for details.
#'   In addition, it contains elements `call` and `optimizer`.
#'
#' @details
#' The `formula` typically looks like:
#'
#' `FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID)`
#'
#' which specifies response and covariates as usual, and exactly one special term
#' defines which covariance structure is used and what are the visit and
#' subject variables.
#'
#' Always use only the first optimizer if multiple optimizers are provided.
#'
#' @export
#'
#' @examples
#' formula <- FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID)
#' data <- fev_data
#' system.time(result <- fit_mmrm(formula, data, rep(1, nrow(fev_data))))
fit_mmrm <- function(formula,
                     data,
                     weights,
                     reml = TRUE,
                     covariance = NULL,
                     tmb_data,
                     formula_parts,
                     control = mmrm_control()) {
  if (missing(formula_parts) || missing(tmb_data)) {
    covariance <- h_reconcile_cov_struct(formula, covariance)
    formula_parts <- h_mmrm_tmb_formula_parts(formula, covariance)

    if (!formula_parts$is_spatial && !is.factor(data[[formula_parts$visit_var]])) {
      stop("Time variable must be a factor for non-spatial covariance structures")
    }

    assert_class(control, "mmrm_control")
    assert_list(control$optimizers, min.len = 1)
    assert_numeric(weights, any.missing = FALSE)
    assert_true(all(weights > 0))
    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
    )
  } else {
    assert_class(tmb_data, "mmrm_tmb_data")
    assert_class(formula_parts, "mmrm_tmb_formula_parts")
  }
  tmb_parameters <- h_mmrm_tmb_parameters(formula_parts, tmb_data, start = control$start, n_groups = tmb_data$n_groups)

  tmb_object <- TMB::MakeADFun(
    data = tmb_data,
    parameters = tmb_parameters,
    hessian = TRUE,
    DLL = "mmrm",
    silent = TRUE
  )
  h_mmrm_tmb_assert_start(tmb_object)
  used_optimizer <- control$optimizers[[1L]]
  used_optimizer_name <- names(control$optimizers)[1L]
  args <- with(
    tmb_object,
    c(
      list(par, fn, gr),
      attr(used_optimizer, "args")
    )
  )
  if (identical(attr(used_optimizer, "use_hessian"), TRUE)) {
    args$hessian <- tmb_object$he
  }
  tmb_opt <- do.call(
    what = used_optimizer,
    args = args
  )
  # Ensure negative log likelihood is stored in `objective` element of list.
  if ("value" %in% names(tmb_opt)) {
    tmb_opt$objective <- tmb_opt$value
    tmb_opt$value <- NULL
  }
  fit <- h_mmrm_tmb_fit(tmb_object, tmb_opt, formula_parts, tmb_data)
  h_mmrm_tmb_check_conv(tmb_opt, fit)
  fit$call <- match.call()
  fit$call$formula <- formula_parts$formula
  fit$optimizer <- used_optimizer_name
  fit
}

#' Obtain List of Jacobian Matrix Entries for Covariance Matrix
#'
#' @description Obtain the Jacobian matrices given the covariance function and variance parameters.
#'
#' @param tmb_data (`mmrm_tmb_data`)\cr produced by [h_mmrm_tmb_data()].
#' @param theta_est (`numeric`)\cr variance parameters point estimate.
#' @param beta_vcov (`matrix`)\cr vairance covariance matrix of coefficients.
#'
#' @return List with one element per variance parameter containing a matrix
#'   of the same dimensions as the covariance matrix. The values are the derivatives
#'   with regards to this variance parameter.
#'
#' @keywords internal
h_jac_list <- function(tmb_data,
                       theta_est,
                       beta_vcov) {
  assert_class(tmb_data, "mmrm_tmb_data")
  assert_numeric(theta_est)
  assert_matrix(beta_vcov)
  .Call(`_mmrm_get_jacobian`, PACKAGE = "mmrm", tmb_data, theta_est, beta_vcov)
}

#' Quadratic Form Calculations
#'
#' @description These helpers are mainly for easier readability and slightly better efficiency
#' of the quadratic forms used in the Satterthwaite calculations.
#'
#' @param center (`matrix`)\cr square numeric matrix with the same dimensions as
#'   `x` as the center of the quadratic form.
#'
#' @name h_quad_form
NULL

#' @describeIn h_quad_form calculates the number `vec %*% center %*% t(vec)`
#'   as a numeric (not a matrix).
#'
#' @param vec (`numeric`)\cr interpreted as a row vector.
#'
#' @keywords internal
h_quad_form_vec <- function(vec, center) {
  vec <- as.vector(vec)
  assert_numeric(vec, any.missing = FALSE)
  assert_matrix(
    center,
    mode = "numeric",
    any.missing = FALSE,
    nrows = length(vec),
    ncols = length(vec)
  )

  sum(vec * (center %*% vec))
}

#' @describeIn h_quad_form calculates the quadratic form `mat %*% center %*% t(mat)`
#'   as a matrix, the result is square and has dimensions identical to the number
#'   of rows in `mat`.
#'
#' @param mat (`matrix`)\cr numeric matrix to be multiplied left and right of
#'   `center`, therefore needs to have as many columns as there are rows and columns
#'   in `center`.
#'
#' @keywords internal
h_quad_form_mat <- function(mat, center) {
  assert_matrix(mat, mode = "numeric", any.missing = FALSE, min.cols = 1L)
  assert_matrix(
    center,
    mode = "numeric",
    any.missing = FALSE,
    nrows = ncol(center),
    ncols = ncol(center)
  )
  mat %*% tcrossprod(center, mat)
}

#' Computation of a Gradient Given Jacobian and Contrast Vector
#'
#' @description Computes the gradient of a linear combination of `beta` given the Jacobian matrix and
#' variance parameters.
#'
#' @param jac_list (`list`)\cr Jacobian list produced e.g. by [h_jac_list()].
#' @param contrast (`numeric`)\cr contrast vector, which needs to have the
#'   same number of elements as there are rows and columns in each element of
#'   `jac_list`.
#'
#' @return Numeric vector which contains the quadratic forms of each element of
#'   `jac_list` with the `contrast` vector.
#'
#' @keywords internal
h_gradient <- function(jac_list, contrast) {
  assert_list(jac_list)
  assert_numeric(contrast)

  vapply(
    jac_list,
    h_quad_form_vec,
    vec = contrast,
    numeric(1L)
  )
}

#' Calculation of Satterthwaite Degrees of Freedom for One-Dimensional Contrast
#'
#' @description Used in [df_1d()] if method is
#' "Satterthwaite".
#'
#' @param object (`mmrm`)\cr the MMRM fit.
#' @param contrast (`numeric`)\cr contrast vector. Note that this should not include
#'   elements for singular coefficient estimates, i.e. only refer to the
#'   actually estimated coefficients.
#'
#' @return List with `est`, `se`, `df`, `t_stat` and `p_val`.
#' @keywords internal
h_df_1d_sat <- function(object, contrast) {
  assert_class(object, "mmrm")
  contrast <- as.numeric(contrast)
  assert_numeric(contrast, len = length(component(object, "beta_est")))

  df <- if (identical(object$vcov, "Asymptotic")) {
    grad <- h_gradient(component(object, "jac_list"), contrast)
    v_num <- 2 * h_quad_form_vec(contrast, component(object, "beta_vcov"))^2
    v_denom <- h_quad_form_vec(grad, component(object, "theta_vcov"))
    v_num / v_denom
  } else if (object$vcov %in% c("Empirical", "Empirical-Jackknife", "Empirical-Bias-Reduced")) {
    contrast_matrix <- Matrix::.bdiag(rep(list(matrix(contrast, nrow = 1)), component(object, "n_subjects")))
    contrast_matrix <- as.matrix(contrast_matrix)
    g_matrix <- h_quad_form_mat(contrast_matrix, object$empirical_df_mat)
    h_tr(g_matrix)^2 / sum(g_matrix^2)
  }

  h_test_1d(object, contrast, df)
}

#' Calculating Denominator Degrees of Freedom for the Multi-Dimensional Case
#'
#' @description Calculates the degrees of freedom for multi-dimensional contrast.
#'
#' @param t_stat_df (`numeric`)\cr `n` t-statistic derived degrees of freedom.
#'
#' @return Usually the calculation is returning `2 * E / (E - n)` where
#'   `E` is the sum of `t / (t - 2)` over all `t_stat_df` values `t`.
#'
#' @note If the input values are two similar to each other then just the average
#'   of them is returned. If any of the inputs is not larger than 2 then 2 is
#'   returned.
#'
#' @keywords internal
h_md_denom_df <- function(t_stat_df) {
  assert_numeric(t_stat_df, min.len = 1L, lower = .Machine$double.xmin, any.missing = FALSE)

  if (test_scalar(t_stat_df)) {
    t_stat_df
  } else if (all(abs(diff(t_stat_df)) < sqrt(.Machine$double.eps))) {
    mean(t_stat_df)
  } else if (any(t_stat_df <= 2)) {
    2
  } else {
    e <- sum(t_stat_df / (t_stat_df - 2))
    2 * e / (e - (length(t_stat_df)))
  }
}

#' Creating F-Statistic Results from One-Dimensional Contrast
#'
#' @description Creates multi-dimensional result from one-dimensional contrast from [df_1d()].
#'
#' @param object (`mmrm`)\cr model fit.
#' @param contrast (`numeric`)\cr one-dimensional contrast.
#'
#' @return The one-dimensional degrees of freedom are calculated and then
#'   based on that the p-value is calculated.
#'
#' @keywords internal
h_df_md_from_1d <- function(object, contrast) {
  res_1d <- h_df_1d_sat(object, contrast)
  list(
    num_df = 1,
    denom_df = res_1d$df,
    f_stat = res_1d$t_stat^2,
    p_val = stats::pf(q = res_1d$t_stat^2, df1 = 1, df2 = res_1d$df, lower.tail = FALSE)
  )
}

#' Calculation of Satterthwaite Degrees of Freedom for Multi-Dimensional Contrast
#'
#' @description Used in [df_md()] if method is "Satterthwaite".
#'
#' @param object (`mmrm`)\cr the MMRM fit.
#' @param contrast (`matrix`)\cr numeric contrast matrix, if given a `numeric`
#'   then this is coerced to a row vector. Note that this should not include
#'   elements for singular coefficient estimates, i.e. only refer to the
#'   actually estimated coefficients.
#'
#' @return List with `num_df`, `denom_df`, `f_stat` and `p_val` (2-sided p-value).
#' @keywords internal
h_df_md_sat <- function(object, contrast) {
  assert_class(object, "mmrm")
  assert_matrix(contrast, mode = "numeric", any.missing = FALSE, ncols = length(component(object, "beta_est")))
  # Early return if we are in the one-dimensional case.
  if (identical(nrow(contrast), 1L)) {
    return(h_df_md_from_1d(object, contrast))
  }

  contrast_cov <- h_quad_form_mat(contrast, component(object, "beta_vcov"))
  eigen_cont_cov <- eigen(contrast_cov)
  eigen_cont_cov_vctrs <- eigen_cont_cov$vectors
  eigen_cont_cov_vals <- eigen_cont_cov$values

  eps <- sqrt(.Machine$double.eps)
  tol <- max(eps * eigen_cont_cov_vals[1], 0)
  rank_cont_cov <- sum(eigen_cont_cov_vals > tol)
  assert_number(rank_cont_cov, lower = .Machine$double.xmin)
  rank_seq <- seq_len(rank_cont_cov)
  vctrs_cont_prod <- crossprod(eigen_cont_cov_vctrs, contrast)[rank_seq, , drop = FALSE]

  # Early return if rank 1.
  if (identical(rank_cont_cov, 1L)) {
    return(h_df_md_from_1d(object, vctrs_cont_prod))
  }

  t_squared_nums <- drop(vctrs_cont_prod %*% object$beta_est)^2
  t_squared_denoms <- eigen_cont_cov_vals[rank_seq]
  t_squared <- t_squared_nums / t_squared_denoms
  f_stat <- sum(t_squared) / rank_cont_cov
  t_stat_df_nums <- 2 * eigen_cont_cov_vals^2
  t_stat_df <- if (identical(object$vcov, "Asymptotic")) {
    grads_vctrs_cont_prod <- lapply(
      rank_seq,
      function(m) h_gradient(component(object, "jac_list"), contrast = vctrs_cont_prod[m, ])
    )
    t_stat_df_denoms <- vapply(
      grads_vctrs_cont_prod,
      h_quad_form_vec,
      center = component(object, "theta_vcov"),
      numeric(1)
    )
    t_stat_df_nums / t_stat_df_denoms
  } else {
    vapply(
      rank_seq,
      function(m) {
        contrast_matrix <- Matrix::.bdiag(
          rep(list(vctrs_cont_prod[m, , drop = FALSE]), component(object, "n_subjects"))
        )
        contrast_matrix <- as.matrix(contrast_matrix)
        g_matrix <- h_quad_form_mat(contrast_matrix, object$empirical_df_mat)
        h_tr(g_matrix)^2 / sum(g_matrix^2)
      },
      FUN.VALUE = 0
    )
  }
  denom_df <- h_md_denom_df(t_stat_df)

  list(
    num_df = rank_cont_cov,
    denom_df = denom_df,
    f_stat = f_stat,
    p_val = stats::pf(q = f_stat, df1 = rank_cont_cov, df2 = denom_df, lower.tail = FALSE)
  )
}

#' Methods for `mmrm_tmb` Objects
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @param object (`mmrm_tmb`)\cr the fitted MMRM object.
#' @param x (`mmrm_tmb`)\cr same as `object`.
#' @param formula (`mmrm_tmb`)\cr same as `object`.
#' @param complete (`flag`)\cr whether to include potential non-estimable
#'   coefficients.
#' @param ... mostly not used;
#'   Exception is `model.matrix()` passing `...` to the default method.
#' @return Depends on the method, see Functions.
#'
#' @name mmrm_tmb_methods
#'
#' @seealso [`mmrm_methods`], [`mmrm_tidiers`] for additional methods.
#'
#' @examples
#' formula <- FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID)
#' object <- fit_mmrm(formula, fev_data, weights = rep(1, nrow(fev_data)))
NULL

#' @describeIn mmrm_tmb_methods obtains the estimated coefficients.
#' @importFrom stats coef
#' @exportS3Method
#' @examples
#' # Estimated coefficients:
#' coef(object)
coef.mmrm_tmb <- function(object, complete = TRUE, ...) {
  assert_flag(complete)
  nm <- if (complete) "beta_est_complete" else "beta_est"
  component(object, name = nm)
}

#' @describeIn mmrm_tmb_methods obtains the fitted values.
#' @importFrom stats fitted
#' @exportS3Method
#' @examples
#' # Fitted values:
#' fitted(object)
fitted.mmrm_tmb <- function(object, ...) {
  fitted_col <- component(object, "x_matrix") %*% component(object, "beta_est")
  fitted_col[, 1L, drop = TRUE]
}

#' @describeIn mmrm_tmb_methods predict conditional means for new data;
#'  optionally with standard errors and confidence or prediction intervals.
#'  Returns a vector of predictions if `se.fit == FALSE` and
#'  `interval == "none"`; otherwise it returns a data.frame with multiple
#'  columns and one row per input data row.
#'
#' @param newdata (`data.frame`)\cr optional new data, otherwise data from `object` is used.
#' @param se.fit (`flag`)\cr indicator if standard errors are required.
#' @param interval (`string`)\cr type of interval calculation. Can be abbreviated.
#' @param level (`number`)\cr tolerance/confidence level.
#' @param nsim (`count`)\cr number of simulations to use.
#' @param conditional (`flag`)\cr indicator if the prediction is conditional on the observation or not.
#'
#' @importFrom stats predict
#' @exportS3Method
#'
#' @examples
#' predict(object, newdata = fev_data)
predict.mmrm_tmb <- function(object,
                             newdata,
                             se.fit = FALSE, # nolint
                             interval = c("none", "confidence", "prediction"),
                             level = 0.95,
                             nsim = 1000L,
                             conditional = FALSE,
                             ...) {
  if (missing(newdata)) {
    newdata <- object$data
  }
  assert_data_frame(newdata)
  orig_row_names <- row.names(newdata)
  assert_flag(se.fit)
  assert_number(level, lower = 0, upper = 1)
  assert_count(nsim, positive = TRUE)
  assert_flag(conditional)
  interval <- match.arg(interval)
  formula_parts <- object$formula_parts
  if (any(object$tmb_data$x_cols_aliased)) {
    warning(
      "In fitted object there are co-linear variables and therefore dropped terms, ",
      "and this could lead to incorrect prediction on new data."
    )
  }
  colnames <- names(Filter(isFALSE, object$tmb_data$x_cols_aliased))
  if (!conditional && interval %in% c("none", "confidence")) {
    # model.matrix always return a complete matrix (no NA allowed)
    x_mat <- stats::model.matrix(object, data = newdata, use_response = FALSE)[, colnames, drop = FALSE]
    x_mat_full <- matrix(
      NA,
      nrow = nrow(newdata), ncol = ncol(x_mat),
      dimnames = list(row.names(newdata), colnames(x_mat))
    )
    x_mat_full[row.names(x_mat), ] <- x_mat
    predictions <- (x_mat_full %*% component(object, "beta_est"))[, 1]
    predictions_raw <- stats::setNames(rep(NA_real_, nrow(newdata)), row.names(newdata))
    predictions_raw[names(predictions)] <- predictions
    if (identical(interval, "none")) {
      return(predictions_raw)
    }
    se <- switch(interval,
      # can be NA if there are aliased cols
      "confidence" = diag(x_mat_full %*% component(object, "beta_vcov") %*% t(x_mat_full)),
      "none" = NA_real_
    )
    res <- cbind(
      fit = predictions, se = se,
      lwr = predictions - stats::qnorm(1 - level / 2) * se, upr = predictions + stats::qnorm(1 - level / 2) * se
    )
    if (!se.fit) {
      res <- res[, setdiff(colnames(res), "se")]
    }
    res_raw <- matrix(
      NA_real_,
      ncol = ncol(res), nrow = nrow(newdata),
      dimnames = list(row.names(newdata), colnames(res))
    )
    res_raw[row.names(res), ] <- res
    return(res_raw)
  }
  tmb_data <- h_mmrm_tmb_data(
    formula_parts, newdata,
    weights = rep(1, nrow(newdata)),
    reml = TRUE,
    singular = "keep",
    drop_visit_levels = FALSE,
    allow_na_response = TRUE,
    drop_levels = FALSE,
    xlev = component(object, "xlev"),
    contrasts = component(object, "contrasts")
  )
  tmb_data$x_matrix <- tmb_data$x_matrix[, colnames, drop = FALSE]
  predictions <- h_get_prediction(
    tmb_data, object$theta_est, object$beta_est, component(object, "beta_vcov")
  )$prediction
  res <- cbind(fit = rep(NA_real_, nrow(newdata)))
  new_order <- match(row.names(tmb_data$full_frame), orig_row_names)
  res[new_order, "fit"] <- predictions[, "fit"]
  se <- switch(interval,
    "confidence" = sqrt(predictions[, "conf_var"]),
    "prediction" = sqrt(h_get_prediction_variance(object, nsim, tmb_data)),
    "none" = NULL
  )
  if (interval != "none") {
    res <- cbind(
      res,
      se = NA_real_
    )
    res[new_order, "se"] <- se
    alpha <- 1 - level
    z <- stats::qnorm(1 - alpha / 2) * res[, "se"]
    res <- cbind(
      res,
      lwr = res[, "fit"] - z,
      upr = res[, "fit"] + z
    )
    if (!se.fit) {
      res <- res[, setdiff(colnames(res), "se")]
    }
  }
  # Use original names.
  row.names(res) <- orig_row_names
  if (ncol(res) == 1) {
    res <- res[, "fit"]
  }
  return(res)
}

#' Get Prediction
#'
#' @description Get predictions with given `data`, `theta`, `beta`, `beta_vcov`.
#'
#' @details See `predict` function in `predict.cpp` which is called internally.
#'
#' @param tmb_data (`mmrm_tmb_data`)\cr object.
#' @param theta (`numeric`)\cr theta value.
#' @param beta (`numeric`)\cr beta value.
#' @param beta_vcov (`matrix`)\cr beta_vcov matrix.
#'
#' @return List with:
#' - `prediction`: Matrix with columns `fit`, `conf_var`, and `var`.
#' - `covariance`: List with subject specific covariance matrices.
#' - `index`: List of zero-based subject indices.
#'
#' @keywords internal
h_get_prediction <- function(tmb_data, theta, beta, beta_vcov) {
  assert_class(tmb_data, "mmrm_tmb_data")
  assert_numeric(theta)
  n_beta <- ncol(tmb_data$x_matrix)
  assert_numeric(beta, finite = TRUE, any.missing = FALSE, len = n_beta)
  assert_matrix(beta_vcov, mode = "numeric", any.missing = FALSE, nrows = n_beta, ncols = n_beta)
  .Call(`_mmrm_predict`, PACKAGE = "mmrm", tmb_data, theta, beta, beta_vcov)
}

#' Get Prediction Variance
#'
#' @description Get prediction variance with given fit, `tmb_data` with the Monte Carlo sampling method.
#'
#' @param object (`mmrm_tmb`)\cr the fitted MMRM.
#' @param nsim (`count`)\cr number of samples.
#' @param tmb_data (`mmrm_tmb_data`)\cr object.
#'
#' @keywords internal
h_get_prediction_variance <- function(object, nsim, tmb_data) {
  assert_class(object, "mmrm_tmb")
  assert_class(tmb_data, "mmrm_tmb_data")
  assert_count(nsim, positive = TRUE)
  theta_chol <- chol(object$theta_vcov)
  n_theta <- length(object$theta_est)
  res <- replicate(nsim, {
    z <- stats::rnorm(n = n_theta)
    theta_sample <- object$theta_est + z %*% theta_chol
    cond_beta_results <- object$tmb_object$report(theta_sample)
    beta_mean <- cond_beta_results$beta
    beta_cov <- cond_beta_results$beta_vcov
    h_get_prediction(tmb_data, theta_sample, beta_mean, beta_cov)$prediction
  })
  mean_of_var <- rowMeans(res[, "var", ])
  var_of_mean <- apply(res[, "fit", ], 1, stats::var)
  mean_of_var + var_of_mean
}

#' @describeIn mmrm_tmb_methods obtains the model frame.
#' @param data (`data.frame`)\cr object in which to construct the frame.
#' @param include (`character`)\cr names of variable types to include.
#'   Must be `NULL` or one or more of `c("subject_var", "visit_var", "group_var", "response_var")`.
#' @param full (`flag`)\cr indicator whether to return full model frame (deprecated).
#' @param na.action (`string`)\cr na action.
#' @importFrom stats model.frame
#' @exportS3Method
#'
#' @details
#' `include` argument controls the variables the returned model frame will include.
#' Possible options are "response_var", "subject_var", "visit_var" and "group_var", representing the
#' response variable, subject variable, visit variable or group variable.
#' `character` values in new data will always be factorized according to the data in the fit
#' to avoid mismatched in levels or issues in `model.matrix`.
#'
#' @examples
#' # Model frame:
#' model.frame(object)
#' model.frame(object, include = "subject_var")
model.frame.mmrm_tmb <- function(formula, data, include = c("subject_var", "visit_var", "group_var", "response_var"),
                                 full, na.action = "na.omit", ...) { # nolint
  # Construct updated formula and data arguments.
  lst_formula_and_data <-
    h_construct_model_frame_inputs(
      formula = formula,
      data = data,
      include = include,
      full = full
    )
  # Only if include is default (full) and also data is missing, and also na.action is na.omit we will
  # use the model frame from the tmb_data.
  include_choice <- c("subject_var", "visit_var", "group_var", "response_var")
  if (missing(data) && setequal(include, include_choice) && identical(h_get_na_action(na.action), stats::na.omit)) {
    ret <- formula$tmb_data$full_frame
    # Remove weights column.
    ret[, "(weights)"] <- NULL
    ret
  } else {
    # Construct data frame to return to users.
    ret <-
      stats::model.frame(
        formula = lst_formula_and_data$formula,
        data = h_get_na_action(na.action)(lst_formula_and_data$data),
        na.action = na.action,
        xlev = stats::.getXlevels(terms(formula), formula$tmb_data$full_frame)
      )
  }
  ret
}


#' Construction of Model Frame Formula and Data Inputs
#'
#' @description
#' Input formulas are converted from mmrm-style to a style compatible
#' with default [stats::model.frame()] and [stats::model.matrix()] methods.
#'
#' The full formula is returned so we can construct, for example, the
#' `model.frame()` including all columns as well as the requested subset.
#' The full set is used to identify rows to include in the reduced model frame.
#'
#' @param formula (`mmrm`)\cr mmrm fit object.
#' @param data optional data frame that will be
#'   passed to `model.frame()` or `model.matrix()`
#' @param include (`character`)\cr names of variable to include
#' @param full (`flag`)\cr indicator whether to return full model frame (deprecated).
#'
#' @return named list with four elements:
#' - `"formula"`: the formula including the columns requested in the `include=` argument.
#' - `"data"`: a data frame including all columns needed in the formula.
#'   full formula are identical
#' @keywords internal
h_construct_model_frame_inputs <- function(formula,
                                           data,
                                           include,
                                           include_choice = c("subject_var", "visit_var", "group_var", "response_var"),
                                           full) {
  if (!missing(full) && identical(full, TRUE)) {
    lifecycle::deprecate_warn("0.3", "model.frame.mmrm_tmb(full)")
    include <- include_choice
  }

  assert_class(formula, classes = "mmrm_tmb")
  assert_subset(include, include_choice)
  if (missing(data)) {
    data <- formula$data
  }
  assert_data_frame(data)

  drop_response <- !"response_var" %in% include
  add_vars <- unlist(formula$formula_parts[include])
  new_formula <- h_add_terms(formula$formula_parts$model_formula, add_vars, drop_response)

  drop_response_full <- !"response_var" %in% include_choice
  add_vars_full <- unlist(formula$formula_parts[include_choice])
  new_formula_full <-
    h_add_terms(formula$formula_parts$model_formula, add_vars_full, drop_response_full)

  # Update data based on the columns in the full formula return.
  all_vars <- all.vars(new_formula_full)
  assert_names(colnames(data), must.include = all_vars)
  data <- data[, all_vars, drop = FALSE]

  # Return list with updated formula, data.
  list(
    formula = new_formula,
    data = data
  )
}

#' @describeIn mmrm_tmb_methods obtains the model matrix.
#' @exportS3Method
#' @param use_response (`flag`)\cr whether to use the response for complete rows.
#'
#' @examples
#' # Model matrix:
#' model.matrix(object)
model.matrix.mmrm_tmb <- function(object, data, use_response = TRUE, ...) { # nolint
  # Always return the utilized model matrix if data not provided.
  if (missing(data)) {
    return(object$tmb_data$x_matrix)
  }
  stats::model.matrix(
    h_add_terms(object$formula_parts$model_formula, NULL, drop_response = !use_response),
    data = data,
    contrasts.arg = attr(object$tmb_data$x_matrix, "contrasts"),
    xlev = component(object, "xlev"),
    ...
  )
}

#' @describeIn mmrm_tmb_methods obtains the terms object.
#' @importFrom stats model.frame
#' @exportS3Method
#'
#' @examples
#' # terms:
#' terms(object)
#' terms(object, include = "subject_var")
terms.mmrm_tmb <- function(x, include = "response_var", ...) { # nolint
  # Construct updated formula and data arguments.
  lst_formula_and_data <-
    h_construct_model_frame_inputs(
      formula = x,
      include = include
    )

  # Use formula method for `terms()` to construct the mmrm terms object.
  stats::terms(
    x = lst_formula_and_data$formula,
    data = lst_formula_and_data$data
  )
}


#' @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, ...) {
  -component(object, "neg_log_lik")
}

#' @describeIn mmrm_tmb_methods obtains the used formula.
#' @importFrom stats formula
#' @exportS3Method
#' @examples
#' # Formula which was used:
#' formula(object)
formula.mmrm_tmb <- function(x, ...) {
  x$formula_parts$formula
}

#' @describeIn mmrm_tmb_methods obtains the variance-covariance matrix estimate
#'   for the coefficients.
#' @importFrom stats vcov
#' @exportS3Method
#' @examples
#' # Variance-covariance matrix estimate for coefficients:
#' vcov(object)
vcov.mmrm_tmb <- function(object, complete = TRUE, ...) {
  assert_flag(complete)
  nm <- if (complete) "beta_vcov_complete" else "beta_vcov"
  component(object, name = nm)
}

#' @describeIn mmrm_tmb_methods obtains the variance-covariance matrix estimate
#'   for the residuals.
#' @param sigma cannot be used (this parameter does not exist in MMRM).
#' @importFrom nlme VarCorr
#' @export VarCorr
#' @aliases VarCorr
#' @exportS3Method
#' @examples
#' # Variance-covariance matrix estimate for residuals:
#' VarCorr(object)
VarCorr.mmrm_tmb <- function(x, sigma = NA, ...) { # nolint
  assert_scalar_na(sigma)

  component(x, name = "varcor")
}

#' @describeIn mmrm_tmb_methods obtains the deviance, which is defined here
#'   as twice the negative log likelihood, which can either be integrated
#'   over the coefficients for REML fits or the usual one for ML fits.
#' @importFrom stats deviance
#' @exportS3Method
#' @examples
#' # REML criterion (twice the negative log likelihood):
#' deviance(object)
deviance.mmrm_tmb <- function(object, ...) {
  2 * component(object, "neg_log_lik")
}

#' @describeIn mmrm_tmb_methods obtains the Akaike Information Criterion,
#'   where the degrees of freedom are the number of variance parameters (`n_theta`).
#'   If `corrected`, then this is multiplied with `m / (m - n_theta - 1)` where
#'   `m` is the number of observations minus the number of coefficients, or
#'   `n_theta + 2` if it is smaller than that \insertCite{hurvich1989regression,burnham1998practical}{mmrm}.
#' @param corrected (`flag`)\cr whether corrected AIC should be calculated.
#' @param k (`number`)\cr the penalty per parameter to be used; default `k = 2`
#'   is the classical AIC.
#' @importFrom stats AIC
#' @exportS3Method
#' @examples
#' # AIC:
#' AIC(object)
#' AIC(object, corrected = TRUE)
#' @references
#' - \insertRef{hurvich1989regression}{mmrm}
#' - \insertRef{burnham1998practical}{mmrm}
AIC.mmrm_tmb <- function(object, corrected = FALSE, ..., k = 2) {
  # nolint
  assert_flag(corrected)
  assert_number(k, lower = 1)

  n_theta <- length(component(object, "theta_est"))
  df <- if (!corrected) {
    n_theta
  } else {
    n_obs <- length(component(object, "y_vector"))
    n_beta <- length(component(object, "beta_est"))
    m <- max(n_theta + 2, n_obs - n_beta)
    n_theta * (m / (m - n_theta - 1))
  }

  2 * component(object, "neg_log_lik") + k * df
}

#' @describeIn mmrm_tmb_methods obtains the Bayesian Information Criterion,
#'   which is using the natural logarithm of the number of subjects for the
#'   penalty parameter `k`.
#' @importFrom stats BIC
#' @exportS3Method
#' @examples
#' # BIC:
#' BIC(object)
BIC.mmrm_tmb <- function(object, ...) {
  # nolint
  k <- log(component(object, "n_subjects"))
  AIC(object, corrected = FALSE, k = k)
}


#' @describeIn mmrm_tmb_methods prints the object.
#' @exportS3Method
print.mmrm_tmb <- function(x,
                           ...) {
  cat("mmrm fit\n\n")

  h_print_call(
    component(x, "call"), component(x, "n_obs"),
    component(x, "n_subjects"), component(x, "n_timepoints")
  )
  h_print_cov(component(x, "cov_type"), component(x, "n_theta"), component(x, "n_groups"))

  cat("Inference:   ")
  cat(ifelse(component(x, "reml"), "REML", "ML"))
  cat("\n")
  cat("Deviance:    ")
  cat(deviance(x))

  cat("\n\nCoefficients: ")
  n_singular_coefs <- sum(component(x, "beta_aliased"))
  if (n_singular_coefs > 0) {
    cat("(", n_singular_coefs, " not defined because of singularities)", sep = "")
  }
  cat("\n")
  print(coef(x, complete = TRUE))

  cat("\nModel Inference Optimization:")

  cat(ifelse(component(x, "convergence") == 0, "\nConverged", "\nFailed to converge"))
  cat(
    " with code", component(x, "convergence"),
    "and message:",
    if (is.null(component(x, "conv_message"))) "No message provided." else tolower(component(x, "conv_message"))
  )
  cat("\n")
  invisible(x)
}


#' @describeIn mmrm_tmb_methods to obtain residuals - either unscaled ('response'), 'pearson' or 'normalized'.
#' @param type (`string`)\cr unscaled (`response`), `pearson` or `normalized`. Default is `response`,
#' and this is the only type available for use with models with a spatial covariance structure.
#' @importFrom stats residuals
#' @exportS3Method
#' @examples
#' # residuals:
#' residuals(object, type = "response")
#' residuals(object, type = "pearson")
#' residuals(object, type = "normalized")
#' @references
#' - \insertRef{galecki2013linear}{mmrm}
residuals.mmrm_tmb <- function(object, type = c("response", "pearson", "normalized"), ...) {
  type <- match.arg(type)
  switch(type,
    "response" = h_residuals_response(object),
    "pearson" = h_residuals_pearson(object),
    "normalized" = h_residuals_normalized(object)
  )
}
#' Calculate Pearson Residuals
#'
#' This is used by [residuals.mmrm_tmb()] to calculate Pearson residuals.
#'
#' @param object (`mmrm_tmb`)\cr the fitted MMRM.
#'
#' @return Vector of residuals.
#'
#' @keywords internal
h_residuals_pearson <- function(object) {
  assert_class(object, "mmrm_tmb")
  h_residuals_response(object) * object$tmb_object$report()$diag_cov_inv_sqrt
}

#' Calculate normalized residuals
#'
#' This is used by [residuals.mmrm_tmb()] to calculate normalized / scaled residuals.
#'
#' @param object (`mmrm_tmb`)\cr the fitted MMRM.
#'
#' @return Vector of residuals
#'
#' @keywords internal
h_residuals_normalized <- function(object) {
  assert_class(object, "mmrm_tmb")
  object$tmb_object$report()$epsilonTilde
}
#' Calculate response residuals.
#'
#' This is used by [residuals.mmrm_tmb()] to calculate response residuals.
#'
#' @param object (`mmrm_tmb`)\cr the fitted MMRM.
#'
#' @return Vector of residuals
#'
#' @keywords internal
h_residuals_response <- function(object) {
  assert_class(object, "mmrm_tmb")
  component(object, "y_vector") - unname(fitted(object))
}

#' @describeIn mmrm_tmb_methods simulate responses from a fitted model according
#'   to the simulation `method`, returning a `data.frame` of dimension `[n, m]`
#'   where n is the number of rows in `newdata`,
#'   and m is the number `nsim` of simulated responses.
#'
#' @param seed unused argument from [stats::simulate()].
#' @param method (`string`)\cr simulation method to use. If "conditional",
#'   simulated values are sampled given the estimated covariance matrix of `object`.
#'   If "marginal", the variance of the estimated covariance matrix is taken into account.
#'
#' @importFrom stats simulate
#' @exportS3Method
simulate.mmrm_tmb <- function(object,
                              nsim = 1,
                              seed = NULL,
                              newdata,
                              ...,
                              method = c("conditional", "marginal")) {
  assert_count(nsim, positive = TRUE)
  assert_null(seed)
  if (missing(newdata)) {
    newdata <- object$data
  }
  assert_data_frame(newdata)
  method <- match.arg(method)


  tmb_data <- h_mmrm_tmb_data(
    object$formula_parts, newdata,
    weights = rep(1, nrow(newdata)),
    reml = TRUE,
    singular = "keep",
    drop_visit_levels = FALSE,
    allow_na_response = TRUE,
    drop_levels = FALSE,
    xlev = component(object, "xlev"),
    contrasts = component(object, "contrasts")
  )
  ret <- if (method == "conditional") {
    predict_res <- h_get_prediction(tmb_data, object$theta_est, object$beta_est, object$beta_vcov)
    as.data.frame(h_get_sim_per_subj(predict_res, tmb_data$n_subjects, nsim))
  } else if (method == "marginal") {
    theta_chol <- t(chol(object$theta_vcov))
    n_theta <- length(object$theta_est)
    as.data.frame(
      sapply(seq_len(nsim), function(x) {
        newtheta <- object$theta_est + theta_chol %*% matrix(stats::rnorm(n_theta), ncol = 1)
        # Recalculate betas with sampled thetas.
        hold <- object$tmb_object$report(newtheta)
        # Resample betas given new beta distribution.
        # We first solve L^\top w = D^{-1/2}z_{sample}:
        w_sample <- backsolve(
          r = hold$XtWX_L,
          x = stats::rnorm(length(hold$beta)) / sqrt(hold$XtWX_D),
          upper.tri = FALSE,
          transpose = TRUE
        )
        # Then we add the mean vector, the beta estimate.
        beta_sample <- hold$beta + w_sample
        predict_res <- h_get_prediction(tmb_data, newtheta, beta_sample, hold$beta_vcov)
        h_get_sim_per_subj(predict_res, tmb_data$n_subjects, 1L)
      })
    )
  }
  orig_row_names <- row.names(newdata)
  new_order <- match(orig_row_names, row.names(tmb_data$full_frame))
  ret[new_order, , drop = FALSE]
}

#' Get simulated values by patient.
#'
#' @param predict_res (`list`)\cr from [h_get_prediction()].
#' @param nsub (`count`)\cr number of subjects.
#' @param nsim (`count`)\cr number of values to simulate.
#'
#' @keywords internal
h_get_sim_per_subj <- function(predict_res, nsub, nsim) {
  assert_list(predict_res)
  assert_count(nsub, positive = TRUE)
  assert_count(nsim, positive = TRUE)

  ret <- matrix(
    predict_res$prediction[, "fit"],
    ncol = nsim,
    nrow = nrow(predict_res$prediction)
  )
  for (i in seq_len(nsub)) {
    # Skip subjects which are not included in predict_res.
    if (length(predict_res$index[[i]]) > 0) {
      # Obtain indices of data.frame belonging to subject i
      # (increment by 1, since indices from cpp are 0-order).
      inds <- predict_res$index[[i]] + 1
      obs <- length(inds)

      # Get relevant covariance matrix for subject i.
      covmat_i <- predict_res$covariance[[i]]
      theta_chol <- t(chol(covmat_i))

      # Simulate epsilon from covariance matrix.
      mus <- ret[inds, , drop = FALSE]
      epsilons <- theta_chol %*% matrix(stats::rnorm(nsim * obs), ncol = nsim)
      ret[inds, ] <- mus + epsilons
    }
  }

  ret
}

#' Calculation of Degrees of Freedom for One-Dimensional Contrast
#'
#' @description `r lifecycle::badge("stable")`
#' Calculates the estimate, adjusted standard error, degrees of freedom,
#' t statistic and p-value for one-dimensional contrast.
#'
#' @param object (`mmrm`)\cr the MMRM fit.
#' @param contrast (`numeric`)\cr contrast vector. Note that this should not include
#'   elements for singular coefficient estimates, i.e. only refer to the
#'   actually estimated coefficients.
#' @return List with `est`, `se`, `df`, `t_stat` and `p_val`.
#' @export
#'
#' @examples
#' object <- mmrm(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data
#' )
#' contrast <- numeric(length(object$beta_est))
#' contrast[3] <- 1
#' df_1d(object, contrast)
df_1d <- function(object, contrast) {
  assert_class(object, "mmrm")
  assert_numeric(contrast, len = length(component(object, "beta_est")), any.missing = FALSE)
  contrast <- as.vector(contrast)
  switch(object$method,
    "Satterthwaite" = h_df_1d_sat(object, contrast),
    "Kenward-Roger" = h_df_1d_kr(object, contrast),
    "Residual" = h_df_1d_res(object, contrast),
    "Between-Within" = h_df_1d_bw(object, contrast),
    stop("Unrecognized degrees of freedom method: ", object$method)
  )
}


#' Calculation of Degrees of Freedom for Multi-Dimensional Contrast
#'
#' @description `r lifecycle::badge("stable")`
#' Calculates the estimate, standard error, degrees of freedom,
#' t statistic and p-value for one-dimensional contrast, depending on the method
#' used in [mmrm()].
#'
#' @param object (`mmrm`)\cr the MMRM fit.
#' @param contrast (`matrix`)\cr numeric contrast matrix, if given a `numeric`
#'   then this is coerced to a row vector. Note that this should not include
#'   elements for singular coefficient estimates, i.e. only refer to the
#'   actually estimated coefficients.
#'
#' @return List with `num_df`, `denom_df`, `f_stat` and `p_val` (2-sided p-value).
#' @export
#'
#' @examples
#' object <- mmrm(
#'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
#'   data = fev_data
#' )
#' contrast <- matrix(data = 0, nrow = 2, ncol = length(object$beta_est))
#' contrast[1, 2] <- contrast[2, 3] <- 1
#' df_md(object, contrast)
df_md <- function(object, contrast) {
  assert_class(object, "mmrm")
  assert_numeric(contrast, any.missing = FALSE)
  if (!is.matrix(contrast)) {
    contrast <- matrix(contrast, ncol = length(contrast))
  }
  assert_matrix(contrast, ncols = length(component(object, "beta_est")))
  if (nrow(contrast) == 0) {
    return(
      list(
        num_df = 0,
        denom_df = NA_real_,
        f_stat = NA_real_,
        p_val = NA_real_
      )
    )
  }
  switch(object$method,
    "Satterthwaite" = h_df_md_sat(object, contrast),
    "Kenward-Roger" = h_df_md_kr(object, contrast),
    "Residual" = h_df_md_res(object, contrast),
    "Between-Within" = h_df_md_bw(object, contrast),
    stop("Unrecognized degrees of freedom method: ", object$method)
  )
}

#' Creating T-Statistic Test Results For One-Dimensional Contrast
#'
#' @description Creates a list of results for one-dimensional contrasts using
#' a t-test statistic and the given degrees of freedom.
#'
#' @inheritParams df_1d
#' @param df (`number`)\cr degrees of freedom for the one-dimensional contrast.
#'
#' @return List with `est`, `se`, `df`, `t_stat` and `p_val` (2-sided p-value).
#'
#' @keywords internal
h_test_1d <- function(object,
                      contrast,
                      df) {
  assert_class(object, "mmrm")
  assert_numeric(contrast, len = length(component(object, "beta_est")))
  assert_number(df, lower = .Machine$double.xmin)

  est <- sum(contrast * component(object, "beta_est"))
  var <- h_quad_form_vec(contrast, component(object, "beta_vcov"))
  se <- sqrt(var)
  t_stat <- est / se
  p_val <- 2 * stats::pt(q = abs(t_stat), df = df, lower.tail = FALSE)

  list(
    est = est,
    se = se,
    df = df,
    t_stat = t_stat,
    p_val = p_val
  )
}

#' Creating F-Statistic Test Results For Multi-Dimensional Contrast
#'
#' @description Creates a list of results for multi-dimensional contrasts using
#' an F-test statistic and the given degrees of freedom.
#'
#' @inheritParams df_md
#' @param contrast (`matrix`)\cr numeric contrast matrix.
#' @param df (`number`)\cr denominator degrees of freedom for the multi-dimensional contrast.
#' @param f_stat_factor (`number`)\cr optional scaling factor on top of the standard F-statistic.
#'
#' @return List with `num_df`, `denom_df`, `f_stat` and `p_val` (2-sided p-value).
#'
#' @keywords internal
h_test_md <- function(object,
                      contrast,
                      df,
                      f_stat_factor = 1) {
  assert_class(object, "mmrm")
  assert_matrix(contrast, ncols = length(component(object, "beta_est")))
  num_df <- nrow(contrast)
  assert_number(df, lower = .Machine$double.xmin)
  assert_number(f_stat_factor, lower = .Machine$double.xmin)

  prec_contrast <- solve(h_quad_form_mat(contrast, component(object, "beta_vcov")))
  contrast_est <- component(object, "beta_est") %*% t(contrast)
  f_statistic <- as.numeric(f_stat_factor / num_df * h_quad_form_mat(contrast_est, prec_contrast))
  p_val <- stats::pf(
    q = f_statistic,
    df1 = num_df,
    df2 = df,
    lower.tail = FALSE
  )

  list(
    num_df = num_df,
    denom_df = df,
    f_stat = f_statistic,
    p_val = p_val
  )
}

#' Methods for `mmrm` Objects
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @param object (`mmrm`)\cr the fitted MMRM including Jacobian and call etc.
#' @param ... not used.
#' @return Depends on the method, see Details and Functions.
#'
#' @details
#' While printing the summary of (`mmrm`)\cr object, the following will be displayed:
#' 1. Formula. The formula used in the model.
#' 2. Data. The data used for analysis, including number of subjects, number of valid observations.
#' 3. Covariance. The covariance structure and number of variance parameters.
#' 4. Method. Restricted maximum likelihood(REML) or maximum likelihood(ML).
#' 5. Model selection criteria. AIC, BIC, log likelihood and deviance.
#' 6. Coefficients. Coefficients of the covariates.
#' 7. Covariance estimate. The covariance estimate(for each group).
#'    1. If the covariance structure is non-spatial, the covariance matrix of all categorical time points available
#'       in data will be displayed.
#'    2. If the covariance structure is spatial, the covariance matrix of two time points with unit distance
#'       will be displayed.
#'
#' `confint` is used to obtain the confidence intervals for the coefficients.
#' Please note that this is different from the confidence interval of difference
#' of least square means from `emmeans`.
#'
#' @name mmrm_methods
#'
#' @seealso [`mmrm_tmb_methods`], [`mmrm_tidiers`] for additional methods.
#'
#' @examples
#' formula <- FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID)
#' object <- mmrm(formula, fev_data)
NULL

#' Coefficients Table for MMRM Fit
#'
#' This is used by [summary.mmrm()] to obtain the coefficients table.
#'
#' @param object (`mmrm`)\cr model fit.
#'
#' @return Matrix with one row per coefficient and columns
#'   `Estimate`, `Std. Error`, `df`, `t value` and `Pr(>|t|)`.
#'
#' @keywords internal
h_coef_table <- function(object) {
  assert_class(object, "mmrm")

  coef_est <- component(object, "beta_est")
  coef_contrasts <- diag(x = rep(1, length(coef_est)))
  rownames(coef_contrasts) <- names(coef_est)
  coef_table <- t(apply(
    coef_contrasts,
    MARGIN = 1L,
    FUN = function(contrast) unlist(df_1d(object, contrast))
  ))
  assert_names(
    colnames(coef_table),
    identical.to = c("est", "se", "df", "t_stat", "p_val")
  )
  colnames(coef_table) <- c("Estimate", "Std. Error", "df", "t value", "Pr(>|t|)")

  coef_aliased <- component(object, "beta_aliased")
  if (any(coef_aliased)) {
    names_coef_na <- names(which(coef_aliased))
    coef_na_table <- matrix(
      data = NA,
      nrow = length(names_coef_na),
      ncol = ncol(coef_table),
      dimnames = list(names_coef_na, colnames(coef_table))
    )
    coef_table <- rbind(coef_table, coef_na_table)[names(coef_aliased), ]
  }

  coef_table
}

#' @describeIn mmrm_methods summarizes the MMRM fit results.
#' @exportS3Method
#' @examples
#' # Summary:
#' summary(object)
summary.mmrm <- function(object, ...) {
  aic_list <- list(
    AIC = AIC(object),
    BIC = BIC(object),
    logLik = logLik(object),
    deviance = deviance(object)
  )
  coefficients <- h_coef_table(object)
  call <- stats::getCall(object)
  components <- component(object, c(
    "cov_type", "reml", "n_groups", "n_theta",
    "n_subjects", "n_timepoints", "n_obs",
    "beta_vcov", "varcor"
  ))
  components$method <- object$method
  components$vcov <- object$vcov
  structure(
    c(
      components,
      list(
        coefficients = coefficients,
        n_singular_coefs = sum(component(object, "beta_aliased")),
        aic_list = aic_list,
        call = call
      )
    ),
    class = "summary.mmrm"
  )
}

#' Printing MMRM Function Call
#'
#' This is used in [print.summary.mmrm()].
#'
#' @param call (`call`)\cr original [mmrm()] function call.
#' @param n_obs (`int`)\cr number of observations.
#' @param n_subjects (`int`)\cr number of subjects.
#' @param n_timepoints (`int`)\cr number of timepoints.
#'
#' @keywords internal
h_print_call <- function(call, n_obs, n_subjects, n_timepoints) {
  pass <- 0
  if (!is.null(tmp <- call$formula)) {
    cat("Formula:    ", deparse(tmp), fill = TRUE)
    rhs <- tmp[[2]]
    pass <- nchar(deparse(rhs))
  }
  if (!is.null(call$data)) {
    cat(
      "Data:       ", deparse(call$data), "(used", n_obs, "observations from",
      n_subjects, "subjects with maximum", n_timepoints, "timepoints)",
      fill = TRUE
    )
  }
  # Display the expression of weights
  if (!is.null(call$weights)) {
    cat("Weights:    ", deparse(call$weights), fill = TRUE)
  }
}

#' Printing MMRM Covariance Type
#'
#' This is used in [print.summary.mmrm()].
#'
#' @param cov_type (`string`)\cr covariance structure abbreviation.
#' @param n_theta (`count`)\cr number of variance parameters.
#' @param n_groups (`count`)\cr number of groups.
#' @keywords internal
h_print_cov <- function(cov_type, n_theta, n_groups) {
  assert_string(cov_type)
  assert_count(n_theta, positive = TRUE)
  assert_count(n_groups, positive = TRUE)
  cov_definition <- switch(cov_type,
    us = "unstructured",
    toep = "Toeplitz",
    toeph = "heterogeneous Toeplitz",
    ar1 = "auto-regressive order one",
    ar1h = "heterogeneous auto-regressive order one",
    ad = "ante-dependence",
    adh = "heterogeneous ante-dependence",
    cs = "compound symmetry",
    csh = "heterogeneous compound symmetry",
    sp_exp = "spatial exponential"
  )

  catstr <- sprintf(
    "Covariance:  %s (%d variance parameters%s)\n",
    cov_definition,
    n_theta,
    ifelse(n_groups == 1L, "", sprintf(" of %d groups", n_groups))
  )
  cat(catstr)
}

#' Printing AIC and other Model Fit Criteria
#'
#' This is used in [print.summary.mmrm()].
#'
#' @param aic_list (`list`)\cr list as part of from [summary.mmrm()].
#' @param digits (`number`)\cr number of decimal places used with [round()].
#'
#' @keywords internal
h_print_aic_list <- function(aic_list,
                             digits = 1) {
  diag_vals <- round(unlist(aic_list), digits)
  diag_vals <- format(diag_vals)
  print(diag_vals, quote = FALSE)
}

#' @describeIn mmrm_methods prints the MMRM fit summary.
#' @exportS3Method
#' @keywords internal
print.summary.mmrm <- function(x,
                               digits = max(3, getOption("digits") - 3),
                               signif.stars = getOption("show.signif.stars"), # nolint
                               ...) {
  cat("mmrm fit\n\n")
  h_print_call(x$call, x$n_obs, x$n_subjects, x$n_timepoints)
  h_print_cov(x$cov_type, x$n_theta, x$n_groups)
  cat("Method:      ", x$method, "\n", sep = "")
  cat("Vcov Method: ", x$vcov, "\n", sep = "")
  cat("Inference:   ")
  cat(ifelse(x$reml, "REML", "ML"))
  cat("\n\n")
  cat("Model selection criteria:\n")
  h_print_aic_list(x$aic_list)
  cat("\n")
  cat("Coefficients: ")
  if (x$n_singular_coefs > 0) {
    cat("(", x$n_singular_coefs, " not defined because of singularities)", sep = "")
  }
  cat("\n")
  stats::printCoefmat(
    x$coefficients,
    zap.ind = 3,
    digits = digits,
    signif.stars = signif.stars
  )
  cat("\n")
  cat("Covariance estimate:\n")
  if (is.list(x$varcor)) {
    for (v in names(x$varcor)) {
      cat(sprintf("Group: %s\n", v))
      print(round(x$varcor[[v]], digits = digits))
    }
  } else {
    print(round(x$varcor, digits = digits))
  }
  cat("\n")
  invisible(x)
}


#' @describeIn mmrm_methods obtain the confidence intervals for the coefficients.
#' @exportS3Method
#' @examples
#' # Confidence Interval:
#' confint(object)
confint.mmrm <- function(object, parm, level = 0.95, ...) {
  cf <- coef(object)
  pnames <- names(cf)
  if (missing(parm)) {
    parm <- pnames
  }
  assert(
    check_subset(parm, pnames),
    check_integerish(parm, lower = 1L, upper = length(cf))
  )
  if (is.numeric(parm)) parm <- pnames[parm]
  assert_number(level, lower = 0, upper = 1)
  a <- (1 - level) / 2
  pct <- paste(format(100 * c(a, 1 - a), trim = TRUE, scientific = FALSE, digits = 3), "%")
  coef_table <- h_coef_table(object)
  df <- coef_table[parm, "df"]
  ses <- coef_table[parm, "Std. Error"]
  fac <- stats::qt(a, df = df)
  ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, pct))
  sefac <- ses * fac
  ci[] <- cf[parm] + c(sefac, -sefac)
  ci
}

#' Capture all Output
#'
#' This function silences all warnings, errors & messages and instead returns a list
#' containing the results (if it didn't error), as well as the warnings, errors
#' and messages and divergence signals as character vectors.
#'
#' @param expr (`expression`)\cr to be executed.
#' @param remove (`list`)\cr optional list with elements `warnings`, `errors`,
#'   `messages` which can be character vectors, which will be removed from the
#'   results if specified.
#' @param divergence (`list`)\cr optional list similar as `remove`, but these
#'   character vectors will be moved to the `divergence` result and signal
#'   that the fit did not converge.
#'
#' @return
#' A list containing
#'
#' - `result`: The object returned by `expr` or `list()` if an error was thrown.
#' - `warnings`: `NULL` or a character vector if warnings were thrown.
#' - `errors`: `NULL` or a string if an error was thrown.
#' - `messages`: `NULL` or a character vector if messages were produced.
#' - `divergence`: `NULL` or a character vector if divergence messages were caught.
#'
#' @keywords internal
h_record_all_output <- function(expr,
                                remove = list(),
                                divergence = list()) {
  # Note: We don't need to and cannot assert `expr` here.
  assert_list(remove, types = "character")
  assert_list(divergence, types = "character")
  env <- new.env()
  result <- withCallingHandlers(
    withRestarts(
      expr,
      muffleStop = function(e) structure(e$message, class = "try-error")
    ),
    message = function(m) {
      msg_without_newline <- gsub(m$message, pattern = "\n$", replacement = "")
      env$message <- c(env$message, msg_without_newline)
      invokeRestart("muffleMessage")
    },
    warning = function(w) {
      env$warning <- c(env$warning, w$message)
      invokeRestart("muffleWarning")
    },
    error = function(e) {
      env$error <- c(env$error, e$message)
      invokeRestart("muffleStop", e)
    }
  )
  list(
    result = result,
    warnings = setdiff(env$warning, c(remove$warnings, divergence$warnings)),
    errors = setdiff(env$error, c(remove$errors, divergence$errors)),
    messages = setdiff(env$message, c(remove$messages, divergence$messages)),
    divergence = c(
      intersect(env$warning, divergence$warnings),
      intersect(env$error, divergence$errors),
      intersect(env$message, divergence$messages)
    )
  )
}

#' Trace of a Matrix
#'
#' @description Obtain the trace of a matrix if the matrix is diagonal, otherwise raise an error.
#'
#' @param x (`matrix`)\cr square matrix input.
#'
#' @return The trace of the square matrix.
#'
#' @keywords internal
h_tr <- function(x) {
  if (nrow(x) != ncol(x)) {
    stop("x must be square matrix")
  }
  sum(Matrix::diag(x))
}

#' Split Control List
#'
#' @description Split the [mmrm_control()] object according to its optimizers and use additional arguments
#' to replace the elements in the original object.
#'
#' @param control (`mmrm_control`)\cr object.
#' @param ... additional parameters to update the `control` object.
#'
#' @return A `list` of `mmrm_control` entries.
#' @keywords internal
h_split_control <- function(control, ...) {
  assert_class(control, "mmrm_control")
  l <- length(control$optimizers)
  lapply(seq_len(l), function(i) {
    ret <- utils::modifyList(control, list(...))
    ret$optimizers <- control$optimizers[i]
    ret
  })
}

#' Obtain Optimizer according to Optimizer String Value
#'
#' @description This function creates optimizer functions with arguments.
#'
#' @param optimizer (`character`)\cr names of built-in optimizers to try, subset
#'   of "L-BFGS-B", "BFGS", "CG" and "nlminb".
#' @param optimizer_fun (`function` or `list` of `function`)\cr alternatively to `optimizer`,
#'   an optimizer function or a list of optimizer functions can be passed directly here.
#' @param optimizer_args (`list`)\cr additional arguments for `optimizer_fun`.
#' @param optimizer_control (`list`)\cr passed to argument `control` in `optimizer_fun`.
#'
#' @details
#' If you want to use only the built-in optimizers:
#' - `optimizer` is a shortcut to create a list of built-in optimizer functions
#'   passed to `optimizer_fun`.
#' - Allowed are "L-BFGS-B", "BFGS", "CG" (using [stats::optim()] with corresponding method)
#'   and "nlminb" (using [stats::nlminb()]).
#' - Other arguments should go into `optimizer_args`.
#'
#' If you want to use your own optimizer function:
#' - Make sure that there are three arguments: parameter (start value), objective function
#'   and gradient function are sequentially in the function arguments.
#' - If there are other named arguments in front of these, make sure they are correctly
#'   specified through `optimizer_args`.
#' - If the hessian can be used, please make sure its argument name is `hessian` and
#'   please add attribute `use_hessian = TRUE` to the function,
#'   using `attr(fun, "use_hessian) <- TRUE`.
#'
#' @return Named `list` of optimizers created by [h_partial_fun_args()].
#'
#' @keywords internal
h_get_optimizers <- function(optimizer = c("L-BFGS-B", "BFGS", "CG", "nlminb"),
                             optimizer_fun = h_optimizer_fun(optimizer),
                             optimizer_args = list(),
                             optimizer_control = list()) {
  if ("automatic" %in% optimizer) {
    lifecycle::deprecate_warn(
      when = "0.2.0",
      what = I("\"automatic\" optimizer"),
      details = "please just omit optimizer argument"
    )
    optimizer_fun <- h_optimizer_fun()
  }
  assert(
    test_function(optimizer_fun),
    test_list(optimizer_fun, types = "function", names = "unique")
  )
  if (is.function(optimizer_fun)) {
    optimizer_fun <- list(custom_optimizer = optimizer_fun)
  }
  lapply(optimizer_fun, function(x) {
    do.call(h_partial_fun_args, c(list(fun = x, control = optimizer_control), optimizer_args))
  })
}

#' Obtain Optimizer Function with Character
#' @description Obtain the optimizer function through the character provided.
#' @param optimizer (`character`)\cr vector of optimizers.
#'
#' @return A (`list`)\cr of optimizer functions generated from [h_partial_fun_args()].
#' @keywords internal
h_optimizer_fun <- function(optimizer = c("L-BFGS-B", "BFGS", "CG", "nlminb")) {
  optimizer <- match.arg(optimizer, several.ok = TRUE)
  lapply(stats::setNames(optimizer, optimizer), function(x) {
    switch(x,
      "L-BFGS-B" = h_partial_fun_args(fun = stats::optim, method = x),
      "BFGS" = h_partial_fun_args(fun = stats::optim, method = x),
      "CG" = h_partial_fun_args(fun = stats::optim, method = x),
      "nlminb" = h_partial_fun_args(fun = stats::nlminb, additional_attr = list(use_hessian = TRUE))
    )
  })
}

#' Create Partial Functions
#' @description Creates partial functions with arguments.
#'
#' @param fun (`function`)\cr to be wrapped.
#' @param ... Additional arguments for `fun`.
#' @param additional_attr (`list`)\cr of additional attributes to apply to the result.
#'
#' @details This function add `args` attribute to the original function,
#' and add an extra class `partial` to the function.
#' `args` is the argument for the function, and elements in `...` will override the existing
#' arguments in attribute `args`. `additional_attr` will override the existing attributes.
#'
#' @return Object with S3 class `"partial"`, a `function` with `args` attribute (and possibly more
#' attributes from `additional_attr`).
#' @keywords internal
h_partial_fun_args <- function(fun, ..., additional_attr = list()) {
  assert_function(fun)
  assert_list(additional_attr, names = "unique")
  a_args <- list(...)
  assert_list(a_args, names = "unique")
  args <- attr(fun, "args")
  if (is.null(args)) {
    args <- list()
  }
  do.call(
    structure,
    args = utils::modifyList(
      list(
        .Data = fun,
        args = utils::modifyList(args, a_args),
        class = c("partial", "function")
      ),
      additional_attr
    )
  )
}

#' Obtain Default Covariance Method
#'
#' @description Obtain the default covariance method depending on
#' the degrees of freedom method used.
#'
#' @param method (`string`)\cr degrees of freedom method.
#'
#' @details The default covariance method is different for different degrees of freedom method.
#' For "Satterthwaite" or "Between-Within", "Asymptotic" is returned.
#' For "Kenward-Roger" only, "Kenward-Roger" is returned.
#' For "Residual" only, "Empirical" is returned.
#'
#' @return String of the default covariance method.
#' @keywords internal
h_get_cov_default <- function(method = c("Satterthwaite", "Kenward-Roger", "Residual", "Between-Within")) {
  assert_string(method)
  method <- match.arg(method)
  switch(method,
    "Residual" = "Empirical",
    "Satterthwaite" = "Asymptotic",
    "Kenward-Roger" = "Kenward-Roger",
    "Between-Within" = "Asymptotic"
  )
}

#' Complete `character` Vector Names From Values
#'
#' @param x (`character` or `list`)\cr value whose names should be completed
#'   from element values.
#'
#' @return A named vector or list.
#'
#' @keywords internal
fill_names <- function(x) {
  n <- names(x)
  is_unnamed <- if (is.null(n)) rep_len(TRUE, length(x)) else n == ""
  names(x)[is_unnamed] <- x[is_unnamed]
  x
}

#' Drop Items from an Indexible
#'
#' Drop elements from an indexible object (`vector`, `list`, etc.).
#'
#' @param x Any object that can be consumed by [seq_along()] and indexed by a
#'   logical vector of the same length.
#' @param n (`integer`)\cr the number of terms to drop.
#'
#' @return A subset of `x`.
#'
#' @keywords internal
drop_elements <- function(x, n) {
  x[seq_along(x) > n]
}

#' Ask for Confirmation on Large Visit Levels
#'
#' @description Ask the user for confirmation if there are too many visit levels
#' for non-spatial covariance structure in interactive sessions.
#'
#' @param x (`numeric`)\cr number of visit levels.
#'
#' @return Logical value `TRUE`.
#' @keywords internal
h_confirm_large_levels <- function(x) {
  assert_count(x)
  allowed_lvls <- x <= getOption("mmrm.max_visits", 100)
  if (allowed_lvls) {
    return(TRUE)
  }
  if (!interactive()) {
    stop("Visit levels too large!", call. = FALSE)
  }
  proceed <- utils::askYesNo(
    paste(
      "Visit levels is possibly too large.",
      "This requires large memory. Are you sure to continue?",
      collapse = " "
    )
  )
  if (!identical(proceed, TRUE)) {
    stop("Visit levels too large!", call. = FALSE)
  }
  return(TRUE)
}

#' Default Value on NULL
#' Return default value when first argument is NULL.
#'
#' @param x Object.
#' @param y Object.
#'
#' @details If `x` is NULL, returns `y`. Otherwise return `x`.
#'
#' @keywords internal
h_default_value <- function(x, y) {
  if (is.null(x)) {
    y
  } else {
    x
  }
}

#' Warn on na.action
#' @keywords internal
h_warn_na_action <- function() {
  if (!identical(getOption("na.action"), "na.omit")) {
    warning("na.action is always set to `na.omit` for `mmrm` fit!")
  }
}

#' Obtain `na.action` as Function
#' @keywords internal
h_get_na_action <- function(na_action) {
  if (is.function(na_action) && identical(methods::formalArgs(na_action), c("object", "..."))) {
    return(na_action)
  }
  if (is.character(na_action) && length(na_action) == 1L) {
    assert_subset(na_action, c("na.omit", "na.exclude", "na.fail", "na.pass", "na.contiguous"))
    return(get(na_action, mode = "function", pos = "package:stats"))
  }
}

#' Validate mmrm Formula
#' @param formula (`formula`)\cr to check.
#'
#' @details In mmrm models, `.` is not allowed as it introduces ambiguity of covariates
#' to be used, so it is not allowed to be in formula.
#'
#' @keywords internal
h_valid_formula <- function(formula) {
  assert_formula(formula)
  if ("." %in% all.vars(formula)) {
    stop("`.` is not allowed in mmrm models!")
  }
}

#' Standard Starting Value
#'
#' @description Obtain standard start values.
#'
#' @param cov_type (`string`)\cr name of the covariance structure.
#' @param n_visits (`int`)\cr number of visits.
#' @param n_groups (`int`)\cr number of groups.
#' @param ... not used.
#'
#' @details
#' `std_start` will try to provide variance parameter from identity matrix.
#' However, for `ar1` and `ar1h` the corresponding values are not ideal because the
#' \eqn{\rho} is usually a positive number thus using 0 as starting value can lead to
#' incorrect optimization result, and we use 0.5 as the initial value of \eqn{\rho}.
#'
#' @return A numeric vector of starting values.
#'
#' @export
std_start <- function(cov_type, n_visits, n_groups, ...) {
  assert_string(cov_type)
  assert_subset(cov_type, cov_types(c("abbr", "habbr")))
  assert_int(n_visits, lower = 1L)
  assert_int(n_groups, lower = 1L)
  start_value <- switch(cov_type,
    us = rep(0, n_visits * (n_visits + 1) / 2),
    toep = rep(0, n_visits),
    toeph = rep(0, 2 * n_visits - 1),
    ar1 = c(0, 0.5),
    ar1h = c(rep(0, n_visits), 0.5),
    ad = rep(0, n_visits),
    adh = rep(0, 2 * n_visits - 1),
    cs = rep(0, 2),
    csh = rep(0, n_visits + 1),
    sp_exp = rep(0, 2)
  )
  rep(start_value, n_groups)
}

#' Empirical Starting Value
#'
#' @description Obtain empirical start value for unstructured covariance
#'
#' @param data (`data.frame`)\cr data used for model fitting.
#' @param model_formula (`formula`)\cr the formula in mmrm model without covariance structure part.
#' @param visit_var (`string`)\cr visit variable.
#' @param subject_var (`string`)\cr subject id variable.
#' @param subject_groups (`factor`)\cr subject group assignment.
#' @param ... not used.
#'
#' @details
#' This `emp_start` only works for unstructured covariance structure.
#' It uses linear regression to first obtain the coefficients and use the residuals
#' to obtain the empirical variance-covariance, and it is then used to obtain the
#' starting values.
#'
#' @note `data` is used instead of `full_frame` because `full_frame` is already
#' transformed if model contains transformations, e.g. `log(FEV1) ~ exp(FEV1_BL)` will
#' drop `FEV1` and `FEV1_BL` but add `log(FEV1)` and `exp(FEV1_BL)` in `full_frame`.
#'
#' @return A numeric vector of starting values.
#'
#' @export
em

1		#' Tidying Methods for `mmrm` Objects
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' These methods tidy the estimates from an `mmrm` object into a
6		#' summary.
7		#'
8		#' @param x (`mmrm`)\cr fitted model.
9		#' @param conf.int (`flag`)\cr if `TRUE` columns for the lower (`conf.low`) and upper bounds
10		#' (`conf.high`) of coefficient estimates are included.
11		#' @param conf.level (`number`)\cr defines the range of the optional confidence internal.
12		#' @param newdata (`data.frame` or `NULL`)\cr optional new data frame.
13		#' @param se_fit (`flag`)\cr whether to return standard errors of fit.
14		#' @param interval (`string`)\cr type of interval calculation.
15		#' @param type.residuals (`string`)\cr passed on to [residuals.mmrm_tmb()].
16		#' @param ... only used by `augment()` to pass arguments to the [predict.mmrm_tmb()] method.
17		#'
18		#' @name mmrm_tidiers
19		#' @aliases mmrm_tidiers
20		#'
21		#' @seealso [`mmrm_methods`], [`mmrm_tmb_methods`] for additional methods.
22		#'
23		#' @examples
24		#' fit <- mmrm(
25		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
26		#' data = fev_data
27		#' )
28		NULL
29
30		#' @describeIn mmrm_tidiers derives tidy `tibble` from an `mmrm` object.
31		#' @exportS3Method
32		#' @examples
33		#' # Applying tidy method to return summary table of covariate estimates.
34		#' fit \|> tidy()
35		#' fit \|> tidy(conf.int = TRUE, conf.level = 0.9)
36		tidy.mmrm <- function(x, # nolint
37		conf.int = FALSE, # nolint
38		conf.level = 0.95, # nolint
39		...) {
40	5x	assert_flag(conf.int)
41	5x	assert_number(conf.level, lower = 0, upper = 1)
42	5x	tbl <- tibble::as_tibble(summary(x)$coefficients, rownames = "term")
43	5x	colnames(tbl) <- c("term", "estimate", "std.error", "df", "statistic", "p.value")
44	5x	coefs <- coef(x)
45	5x	if (length(coefs) != nrow(tbl)) {
46	!	coefs <- tibble::enframe(coefs, name = "term", value = "estimate")
47	!	tbl <- merge(coefs, tbl, by = c("term", "estimate"))
48		}
49	5x	if (conf.int) {
50	4x	ci <- h_tbl_confint_terms(x, level = conf.level)
51	4x	tbl <- tibble::as_tibble(merge(tbl, ci, by = "term"))
52		}
53	5x	tbl
54		}
55
56		#' @describeIn mmrm_tidiers derives `glance` `tibble` from an `mmrm` object.
57		#' @exportS3Method
58		#' @examples
59		#' # Applying glance method to return summary table of goodness of fit statistics.
60		#' fit \|> glance()
61		glance.mmrm <- function(x, ...) { # nolint
62	1x	tibble::as_tibble(summary(x)$aic_list)
63		}
64
65		#' @describeIn mmrm_tidiers derives `augment` `tibble` from an `mmrm` object.
66		#' @exportS3Method
67		#' @examples
68		#' # Applying augment method to return merged `tibble` of model data, fitted and residuals.
69		#' fit \|> augment()
70		#' fit \|> augment(interval = "confidence")
71		#' fit \|> augment(type.residuals = "pearson")
72		augment.mmrm <- function(x, # nolint
73		newdata = NULL,
74		interval = c("none", "confidence", "prediction"),
75		se_fit = (interval != "none"),
76		type.residuals = c("response", "pearson", "normalized"), # nolint
77		...) {
78	9x	type.residuals <- match.arg(type.residuals) # nolint
79	9x	resid_df <- NULL
80	9x	if (is.null(newdata)) {
81	4x	newdata <- stats::get_all_vars(x, data = stats::na.omit(x$data))
82	4x	resid_df <- data.frame(
83	4x	.rownames = rownames(newdata),
84	4x	.resid = unname(residuals(x, type = type.residuals))
85		)
86		}
87	9x	interval <- match.arg(interval)
88
89	9x	tbl <- h_newdata_add_pred(
90	9x	x,
91	9x	newdata = newdata,
92	9x	se_fit = se_fit,
93	9x	interval = interval,
94		...
95		)
96	9x	if (!is.null(resid_df)) {
97	4x	tbl <- merge(tbl, resid_df, by = ".rownames")
98	4x	tbl$.rownames <- as.numeric(tbl$.rownames)
99	4x	tbl <- tbl[order(tbl$.rownames), , drop = FALSE]
100		}
101	9x	tibble::as_tibble(tbl)
102		}
103
104		#' Extract `tibble` with Confidence Intervals and Term Names
105		#'
106		#' This is used in [tidy.mmrm()].
107		#'
108		#' @param x (`mmrm`)\cr fit object.
109		#' @param ... passed to [stats::confint()], hence not used at the moment.
110		#'
111		#' @return A `tibble` with `term`, `conf.low`, `conf.high` columns.
112		#'
113		#' @keywords internal
114		h_tbl_confint_terms <- function(x, ...) {
115	8x	df <- stats::confint(x, ...)
116	8x	tbl <- tibble::as_tibble(df, rownames = "term", .name_repair = "minimal")
117	8x	names(tbl) <- c("term", "conf.low", "conf.high")
118	8x	tbl
119		}
120
121		#' Add Prediction Results to New Data
122		#'
123		#' This is used in [augment.mmrm()].
124		#'
125		#' @param x (`mmrm`)\cr fit.
126		#' @param newdata (`data.frame`)\cr data to predict.
127		#' @param se_fit (`flag`)\cr whether to return standard error of prediction,
128		#' can only be used when `interval` is not "none".
129		#' @param interval (`string`)\cr type of interval.
130		#' @param ... passed to [predict.mmrm_tmb()].
131		#'
132		#' @return The `newdata` as a `tibble` with additional columns `.fitted`,
133		#' `.lower`, `.upper` (if interval is not `none`) and `.se.fit` (if `se_fit`
134		#' requested).
135		#'
136		#' @keywords internal
137		h_newdata_add_pred <- function(x,
138		newdata,
139		se_fit,
140		interval,
141		...) {
142	13x	assert_class(x, "mmrm")
143	13x	assert_data_frame(newdata)
144	13x	assert_flag(se_fit)
145	13x	assert_string(interval)
146	13x	if (interval == "none") {
147	7x	assert_false(se_fit)
148		}
149
150	12x	tbl <- h_df_to_tibble(newdata)
151	12x	pred_results <- predict(
152	12x	x,
153	12x	newdata = newdata,
154	12x	na.action = stats::na.pass,
155	12x	se.fit = se_fit,
156	12x	interval = interval,
157		...
158		)
159	12x	if (interval == "none") {
160	6x	assert_numeric(pred_results)
161	6x	tbl$.fitted <- unname(pred_results)
162		} else {
163	6x	assert_matrix(pred_results)
164	6x	tbl$.fitted <- unname(pred_results[, "fit"])
165	6x	tbl$.lower <- unname(pred_results[, "lwr"])
166	6x	tbl$.upper <- unname(pred_results[, "upr"])
167		}
168	12x	if (se_fit) {
169	5x	tbl$.se.fit <- unname(pred_results[, "se"])
170		}
171	12x	tbl
172		}
173
174		#' Coerce a Data Frame to a `tibble`
175		#'
176		#' This is used in [h_newdata_add_pred()].
177		#'
178		#' @details This is only a thin wrapper around [tibble::as_tibble()], except
179		#' giving a useful error message and it checks for `rownames` and adds them
180		#' as a new column `.rownames` if they are not just a numeric sequence as
181		#' per the [tibble::has_rownames()] decision.
182		#'
183		#' @param data (`data.frame`)\cr what to coerce.
184		#'
185		#' @return The `data` as a `tibble`, potentially with a `.rownames` column.
186		#'
187		#' @keywords internal
188		h_df_to_tibble <- function(data) {
189	15x	tryCatch(tbl <- tibble::as_tibble(data), error = function(cnd) {
190	1x	stop("Could not coerce data to `tibble`. Try explicitly passing a",
191	1x	"dataset to either the `data` or `newdata` argument.",
192	1x	call. = FALSE
193		)
194		})
195	14x	if (tibble::has_rownames(data)) {
196	5x	tbl <- tibble::add_column(tbl, .rownames = rownames(data), .before = TRUE)
197		}
198	14x	tbl
199		}

1		#' Calculation of Degrees of Freedom for One-Dimensional Contrast
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#' Calculates the estimate, adjusted standard error, degrees of freedom,
5		#' t statistic and p-value for one-dimensional contrast.
6		#'
7		#' @param object (`mmrm`)\cr the MMRM fit.
8		#' @param contrast (`numeric`)\cr contrast vector. Note that this should not include
9		#' elements for singular coefficient estimates, i.e. only refer to the
10		#' actually estimated coefficients.
11		#' @return List with `est`, `se`, `df`, `t_stat` and `p_val`.
12		#' @export
13		#'
14		#' @examples
15		#' object <- mmrm(
16		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
17		#' data = fev_data
18		#' )
19		#' contrast <- numeric(length(object$beta_est))
20		#' contrast[3] <- 1
21		#' df_1d(object, contrast)
22		df_1d <- function(object, contrast) {
23	338x	assert_class(object, "mmrm")
24	338x	assert_numeric(contrast, len = length(component(object, "beta_est")), any.missing = FALSE)
25	338x	contrast <- as.vector(contrast)
26	338x	switch(object$method,
27	318x	"Satterthwaite" = h_df_1d_sat(object, contrast),
28	19x	"Kenward-Roger" = h_df_1d_kr(object, contrast),
29	!	"Residual" = h_df_1d_res(object, contrast),
30	1x	"Between-Within" = h_df_1d_bw(object, contrast),
31	!	stop("Unrecognized degrees of freedom method: ", object$method)
32		)
33		}
34
35
36		#' Calculation of Degrees of Freedom for Multi-Dimensional Contrast
37		#'
38		#' @description `r lifecycle::badge("stable")`
39		#' Calculates the estimate, standard error, degrees of freedom,
40		#' t statistic and p-value for one-dimensional contrast, depending on the method
41		#' used in [mmrm()].
42		#'
43		#' @param object (`mmrm`)\cr the MMRM fit.
44		#' @param contrast (`matrix`)\cr numeric contrast matrix, if given a `numeric`
45		#' then this is coerced to a row vector. Note that this should not include
46		#' elements for singular coefficient estimates, i.e. only refer to the
47		#' actually estimated coefficients.
48		#'
49		#' @return List with `num_df`, `denom_df`, `f_stat` and `p_val` (2-sided p-value).
50		#' @export
51		#'
52		#' @examples
53		#' object <- mmrm(
54		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
55		#' data = fev_data
56		#' )
57		#' contrast <- matrix(data = 0, nrow = 2, ncol = length(object$beta_est))
58		#' contrast[1, 2] <- contrast[2, 3] <- 1
59		#' df_md(object, contrast)
60		df_md <- function(object, contrast) {
61	150x	assert_class(object, "mmrm")
62	150x	assert_numeric(contrast, any.missing = FALSE)
63	150x	if (!is.matrix(contrast)) {
64	113x	contrast <- matrix(contrast, ncol = length(contrast))
65		}
66	150x	assert_matrix(contrast, ncols = length(component(object, "beta_est")))
67	150x	if (nrow(contrast) == 0) {
68	1x	return(
69	1x	list(
70	1x	num_df = 0,
71	1x	denom_df = NA_real_,
72	1x	f_stat = NA_real_,
73	1x	p_val = NA_real_
74		)
75		)
76		}
77	149x	switch(object$method,
78	145x	"Satterthwaite" = h_df_md_sat(object, contrast),
79	3x	"Kenward-Roger" = h_df_md_kr(object, contrast),
80	!	"Residual" = h_df_md_res(object, contrast),
81	1x	"Between-Within" = h_df_md_bw(object, contrast),
82	!	stop("Unrecognized degrees of freedom method: ", object$method)
83		)
84		}
85
86		#' Creating T-Statistic Test Results For One-Dimensional Contrast
87		#'
88		#' @description Creates a list of results for one-dimensional contrasts using
89		#' a t-test statistic and the given degrees of freedom.
90		#'
91		#' @inheritParams df_1d
92		#' @param df (`number`)\cr degrees of freedom for the one-dimensional contrast.
93		#'
94		#' @return List with `est`, `se`, `df`, `t_stat` and `p_val` (2-sided p-value).
95		#'
96		#' @keywords internal
97		h_test_1d <- function(object,
98		contrast,
99		df) {
100	486x	assert_class(object, "mmrm")
101	486x	assert_numeric(contrast, len = length(component(object, "beta_est")))
102	486x	assert_number(df, lower = .Machine$double.xmin)
103
104	486x	est <- sum(contrast * component(object, "beta_est"))
105	486x	var <- h_quad_form_vec(contrast, component(object, "beta_vcov"))
106	486x	se <- sqrt(var)
107	486x	t_stat <- est / se
108	486x	p_val <- 2 * stats::pt(q = abs(t_stat), df = df, lower.tail = FALSE)
109
110	486x	list(
111	486x	est = est,
112	486x	se = se,
113	486x	df = df,
114	486x	t_stat = t_stat,
115	486x	p_val = p_val
116		)
117		}
118
119		#' Creating F-Statistic Test Results For Multi-Dimensional Contrast
120		#'
121		#' @description Creates a list of results for multi-dimensional contrasts using
122		#' an F-test statistic and the given degrees of freedom.
123		#'
124		#' @inheritParams df_md
125		#' @param contrast (`matrix`)\cr numeric contrast matrix.
126		#' @param df (`number`)\cr denominator degrees of freedom for the multi-dimensional contrast.
127		#' @param f_stat_factor (`number`)\cr optional scaling factor on top of the standard F-statistic.
128		#'
129		#' @return List with `num_df`, `denom_df`, `f_stat` and `p_val` (2-sided p-value).
130		#'
131		#' @keywords internal
132		h_test_md <- function(object,
133		contrast,
134		df,
135		f_stat_factor = 1) {
136	15x	assert_class(object, "mmrm")
137	15x	assert_matrix(contrast, ncols = length(component(object, "beta_est")))
138	15x	num_df <- nrow(contrast)
139	15x	assert_number(df, lower = .Machine$double.xmin)
140	15x	assert_number(f_stat_factor, lower = .Machine$double.xmin)
141
142	15x	prec_contrast <- solve(h_quad_form_mat(contrast, component(object, "beta_vcov")))
143	15x	contrast_est <- component(object, "beta_est") %*% t(contrast)
144	15x	f_statistic <- as.numeric(f_stat_factor / num_df * h_quad_form_mat(contrast_est, prec_contrast))
145	15x	p_val <- stats::pf(
146	15x	q = f_statistic,
147	15x	df1 = num_df,
148	15x	df2 = df,
149	15x	lower.tail = FALSE
150		)
151
152	15x	list(
153	15x	num_df = num_df,
154	15x	denom_df = df,
155	15x	f_stat = f_statistic,
156	15x	p_val = p_val
157		)
158		}

1		#' Fitting an MMRM with Single Optimizer
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' This function helps to fit an MMRM using `TMB` with a single optimizer,
6		#' while capturing messages and warnings.
7		#'
8		#' @inheritParams mmrm
9		#' @param control (`mmrm_control`)\cr object.
10		#' @param tmb_data (`mmrm_tmb_data`)\cr object.
11		#' @param formula_parts (`mmrm_tmb_formula_parts`)\cr object.
12		#' @param ... Additional arguments to pass to [mmrm_control()].
13		#'
14		#' @details
15		#' `fit_single_optimizer` will fit the `mmrm` model using the `control` provided.
16		#' If there are multiple optimizers provided in `control`, only the first optimizer
17		#' will be used.
18		#' If `tmb_data` and `formula_parts` are both provided, `formula`, `data`, `weights`,
19		#' `reml`, and `covariance` are ignored.
20		#'
21		#' @return The `mmrm_fit` object, with additional attributes containing warnings,
22		#' messages, optimizer used and convergence status in addition to the
23		#' `mmrm_tmb` contents.
24		#' @export
25		#'
26		#' @examples
27		#' mod_fit <- fit_single_optimizer(
28		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
29		#' data = fev_data,
30		#' weights = rep(1, nrow(fev_data)),
31		#' optimizer = "nlminb"
32		#' )
33		#' attr(mod_fit, "converged")
34		fit_single_optimizer <- function(formula,
35		data,
36		weights,
37		reml = TRUE,
38		covariance = NULL,
39		tmb_data,
40		formula_parts,
41		...,
42		control = mmrm_control(...)) {
43	199x	to_remove <- list(
44		# Transient visit to invalid parameters.
45	199x	warnings = c("NA/NaN function evaluation")
46		)
47	199x	as_diverged <- list(
48	199x	errors = c(
49	199x	"NA/NaN Hessian evaluation",
50	199x	"L-BFGS-B needs finite values of 'fn'"
51		)
52		)
53	199x	if (missing(tmb_data) \|\| missing(formula_parts)) {
54	14x	h_valid_formula(formula)
55	13x	assert_data_frame(data)
56	13x	assert_numeric(weights, any.missing = FALSE, lower = .Machine$double.xmin)
57	13x	assert_flag(reml)
58	13x	assert_class(control, "mmrm_control")
59	13x	assert_list(control$optimizers, names = "unique", types = c("function", "partial"))
60	13x	quiet_fit <- h_record_all_output(
61	13x	fit_mmrm(
62	13x	formula = formula,
63	13x	data = data,
64	13x	weights = weights,
65	13x	reml = reml,
66	13x	covariance = covariance,
67	13x	control = control
68		),
69	13x	remove = to_remove,
70	13x	divergence = as_diverged
71		)
72		} else {
73	185x	assert_class(tmb_data, "mmrm_tmb_data")
74	185x	assert_class(formula_parts, "mmrm_tmb_formula_parts")
75	185x	quiet_fit <- h_record_all_output(
76	185x	fit_mmrm(
77	185x	formula_parts = formula_parts,
78	185x	tmb_data = tmb_data,
79	185x	control = control
80		),
81	185x	remove = to_remove,
82	185x	divergence = as_diverged
83		)
84		}
85	198x	if (length(quiet_fit$errors)) {
86	4x	stop(quiet_fit$errors)
87		}
88	194x	converged <- (length(quiet_fit$warnings) == 0L) &&
89	194x	(length(quiet_fit$divergence) == 0L) &&
90	194x	isTRUE(quiet_fit$result$opt_details$convergence == 0)
91	194x	structure(
92	194x	quiet_fit$result,
93	194x	warnings = quiet_fit$warnings,
94	194x	messages = quiet_fit$messages,
95	194x	divergence = quiet_fit$divergence,
96	194x	converged = converged,
97	194x	class = c("mmrm_fit", class(quiet_fit$result))
98		)
99		}
100
101		#' Summarizing List of Fits
102		#'
103		#' @param all_fits (`list` of `mmrm_fit` or `try-error`)\cr list of fits.
104		#'
105		#' @return List with `warnings`, `messages`, `log_liks` and `converged` results.
106		#' @keywords internal
107		h_summarize_all_fits <- function(all_fits) {
108	8x	assert_list(all_fits, types = c("mmrm_fit", "try-error"))
109	8x	is_error <- vapply(all_fits, is, logical(1), class2 = "try-error")
110
111	8x	warnings <- messages <- vector(mode = "list", length = length(all_fits))
112	8x	warnings[is_error] <- lapply(all_fits[is_error], as.character)
113	8x	warnings[!is_error] <- lapply(all_fits[!is_error], attr, which = "warnings")
114	8x	messages[!is_error] <- lapply(all_fits[!is_error], attr, which = "messages")
115	8x	log_liks <- as.numeric(rep(NA, length.out = length(all_fits)))
116	8x	log_liks[!is_error] <- vapply(all_fits[!is_error], stats::logLik, numeric(1L))
117	8x	converged <- rep(FALSE, length.out = length(all_fits))
118	8x	converged[!is_error] <- vapply(all_fits[!is_error], attr, logical(1), which = "converged")
119
120	8x	list(
121	8x	warnings = warnings,
122	8x	messages = messages,
123	8x	log_liks = log_liks,
124	8x	converged = converged
125		)
126		}
127
128		#' Refitting MMRM with Multiple Optimizers
129		#'
130		#' @description `r lifecycle::badge("stable")`
131		#'
132		#' @param fit (`mmrm_fit`)\cr original model fit from [fit_single_optimizer()].
133		#' @param ... Additional arguments passed to [mmrm_control()].
134		#' @param control (`mmrm_control`)\cr object.
135		#'
136		#' @return The best (in terms of log likelihood) fit which converged.
137		#'
138		#' @note For Windows, no parallel computations are currently implemented.
139		#' @export
140		#'
141		#' @examples
142		#' fit <- fit_single_optimizer(
143		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
144		#' data = fev_data,
145		#' weights = rep(1, nrow(fev_data)),
146		#' optimizer = "nlminb"
147		#' )
148		#' best_fit <- refit_multiple_optimizers(fit)
149		refit_multiple_optimizers <- function(fit,
150		...,
151		control = mmrm_control(...)) {
152	6x	assert_class(fit, "mmrm_fit")
153	6x	assert_class(control, "mmrm_control")
154
155	6x	n_cores_used <- ifelse(
156	6x	.Platform$OS.type == "windows",
157	6x	1L,
158	6x	min(
159	6x	length(control$optimizers),
160	6x	control$n_cores
161		)
162		)
163	6x	controls <- h_split_control(
164	6x	control,
165	6x	start = fit$theta_est
166		)
167
168		# Take the results from old fit as starting values for new fits.
169	6x	all_fits <- suppressWarnings(parallel::mcmapply(
170	6x	FUN = fit_single_optimizer,
171	6x	control = controls,
172	6x	MoreArgs = list(
173	6x	tmb_data = fit$tmb_data,
174	6x	formula_parts = fit$formula_parts
175		),
176	6x	mc.cores = n_cores_used,
177	6x	mc.silent = TRUE,
178	6x	SIMPLIFY = FALSE
179		))
180	6x	all_fits <- c(all_fits, list(old_result = fit))
181
182		# Find the results that are ok and return best in terms of log-likelihood.
183	6x	all_fits_summary <- h_summarize_all_fits(all_fits)
184	6x	is_ok <- all_fits_summary$converged
185	6x	if (!any(is_ok)) {
186	1x	stop(
187	1x	"No optimizer led to a successful model fit. ",
188	1x	"Please try to use a different covariance structure or other covariates."
189		)
190		}
191	5x	best_optimizer <- which.max(all_fits_summary$log_liks[is_ok])
192	5x	all_fits[[which(is_ok)[best_optimizer]]]
193		}
194
195		#' Control Parameters for Fitting an MMRM
196		#'
197		#' @description `r lifecycle::badge("stable")`
198		#' Fine-grained specification of the MMRM fit details is possible using this
199		#' control function.
200		#'
201		#' @param n_cores (`count`)\cr number of cores to be used.
202		#' @param method (`string`)\cr adjustment method for degrees of freedom.
203		#' @param vcov (`string`)\cr coefficients covariance matrix adjustment method.
204		#' @param start (`NULL`, `numeric` or `function`)\cr optional start values for variance
205		#' parameters. See details for more information.
206		#' @param accept_singular (`flag`)\cr whether singular design matrices are reduced
207		#' to full rank automatically and additional coefficient estimates will be missing.
208		#' @param optimizers (`list`)\cr optimizer specification, created with [h_get_optimizers()].
209		#' @param drop_visit_levels (`flag`)\cr whether to drop levels for visit variable,
210		#' if visit variable is a factor, see details.
211		#' @param ... additional arguments passed to [h_get_optimizers()].
212		#'
213		#' @details
214		# - The `drop_visit_levels` flag will decide whether unobserved visits will be kept for analysis.
215		#' For example, if the data only has observations at visits `VIS1`, `VIS3` and `VIS4`, by default
216		#' they are treated to be equally spaced, the distance from `VIS1` to `VIS3`, and from `VIS3` to `VIS4`,
217		#' are identical. However, you can manually convert this visit into a factor, with
218		#' `levels = c("VIS1", "VIS2", "VIS3", "VIS4")`, and also use `drop_visits_levels = FALSE`,
219		#' then the distance from `VIS1` to `VIS3` will be double, as `VIS2` is a valid visit.
220		#' However, please be cautious because this can lead to convergence failure
221		#' when using an unstructured covariance matrix and there are no observations
222		#' at the missing visits.
223		#' - The `method` and `vcov` arguments specify the degrees of freedom and coefficients
224		#' covariance matrix adjustment methods, respectively.
225		#' - Allowed `vcov` includes: "Asymptotic", "Kenward-Roger", "Kenward-Roger-Linear", "Empirical" (CR0),
226		#' "Empirical-Jackknife" (CR3), and "Empirical-Bias-Reduced" (CR2).
227		#' - Allowed `method` includes: "Satterthwaite", "Kenward-Roger", "Between-Within" and "Residual".
228		#' - If `method` is "Kenward-Roger" then only "Kenward-Roger" or "Kenward-Roger-Linear" are allowed for `vcov`.
229		#' - The `vcov` argument can be `NULL` to use the default covariance method depending on the `method`
230		#' used for degrees of freedom, see the following table:
231		#'
232		#' \| `method` \| Default `vcov`\|
233		#' \|-----------\|----------\|
234		#' \|Satterthwaite\| Asymptotic\|
235		#' \|Kenward-Roger\| Kenward-Roger\|
236		#' \|Residual\| Empirical\|
237		#' \|Between-Within\| Asymptotic\|
238		#'
239		#' - Please note that "Kenward-Roger" for "Unstructured" covariance gives different results
240		#' compared to SAS; Use "Kenward-Roger-Linear" for `vcov` instead for better matching
241		#' of the SAS results.
242		#'
243		#' - The argument `start` is used to facilitate the choice of initial values for fitting the model.
244		#' If `function` is provided, make sure its parameter is a valid element of `mmrm_tmb_data`
245		#' or `mmrm_tmb_formula_parts` and it returns a numeric vector.
246		#' By default or if `NULL` is provided, `std_start` will be used.
247		#' Other implemented methods include `emp_start`.
248		#'
249		#' @return List of class `mmrm_control` with the control parameters.
250		#' @export
251		#'
252		#' @examples
253		#' mmrm_control(
254		#' optimizer_fun = stats::optim,
255		#' optimizer_args = list(method = "L-BFGS-B")
256		#' )
257		mmrm_control <- function(n_cores = 1L,
258		method = c("Satterthwaite", "Kenward-Roger", "Residual", "Between-Within"),
259		vcov = NULL,
260		start = std_start,
261		accept_singular = TRUE,
262		drop_visit_levels = TRUE,
263		...,
264		optimizers = h_get_optimizers(...)) {
265	243x	assert_count(n_cores, positive = TRUE)
266	243x	assert_character(method)
267	243x	if (is.null(start)) {
268	1x	start <- std_start
269		}
270	243x	assert(
271	243x	check_function(start, args = "..."),
272	243x	check_numeric(start, null.ok = FALSE),
273	243x	combine = "or"
274		)
275	243x	assert_flag(accept_singular)
276	243x	assert_flag(drop_visit_levels)
277	243x	assert_list(optimizers, names = "unique", types = c("function", "partial"))
278	243x	assert_string(vcov, null.ok = TRUE)
279	243x	method <- match.arg(method)
280	243x	if (is.null(vcov)) {
281	192x	vcov <- h_get_cov_default(method)
282		}
283	243x	assert_subset(
284	243x	vcov,
285	243x	c(
286	243x	"Asymptotic",
287	243x	"Empirical",
288	243x	"Empirical-Bias-Reduced",
289	243x	"Empirical-Jackknife",
290	243x	"Kenward-Roger",
291	243x	"Kenward-Roger-Linear"
292		)
293		)
294	243x	if (xor(identical(method, "Kenward-Roger"), vcov %in% c("Kenward-Roger", "Kenward-Roger-Linear"))) {
295	5x	stop(paste(
296	5x	"Kenward-Roger degrees of freedom must work together with Kenward-Roger",
297	5x	"or Kenward-Roger-Linear covariance!"
298		))
299		}
300	238x	structure(
301	238x	list(
302	238x	optimizers = optimizers,
303	238x	start = start,
304	238x	accept_singular = accept_singular,
305	238x	method = method,
306	238x	vcov = vcov,
307	238x	n_cores = as.integer(n_cores),
308	238x	drop_visit_levels = drop_visit_levels
309		),
310	238x	class = "mmrm_control"
311		)
312		}
313
314		#' Fit an MMRM
315		#'
316		#' @description `r lifecycle::badge("stable")`
317		#'
318		#' This is the main function fitting the MMRM.
319		#'
320		#' @param formula (`formula`)\cr the model formula, see details.
321		#' @param data (`data`)\cr the data to be used for the model.
322		#' @param weights (`vector`)\cr an optional vector of weights to be used in
323		#' the fitting process. Should be `NULL` or a numeric vector.
324		#' @param reml (`flag`)\cr whether restricted maximum likelihood (REML)
325		#' estimation is used, otherwise maximum likelihood (ML) is used.
326		#' @param covariance (`cov_struct`)\cr a covariance structure type definition
327		#' as produced with [cov_struct()], or value that can be coerced to a
328		#' covariance structure using [as.cov_struct()]. If no value is provided,
329		#' a structure is derived from the provided formula.
330		#' @param control (`mmrm_control`)\cr fine-grained fitting specifications list
331		#' created with [mmrm_control()].
332		#' @param ... arguments passed to [mmrm_control()].
333		#'
334		#' @details
335		#' The `formula` typically looks like:
336		#' `FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID)`
337		#' so specifies response and covariates as usual, and exactly one special term
338		#' defines which covariance structure is used and what are the time point and
339		#' subject variables. The covariance structures in the formula can be
340		#' found in [`covariance_types`].
341		#'
342		#' The time points have to be unique for each subject. That is,
343		#' there cannot be time points with multiple observations for any subject.
344		#' The rationale is that these observations would need to be correlated, but it
345		#' is not possible within the currently implemented covariance structure framework
346		#' to do that correctly. Moreover, for non-spatial covariance structures, the time
347		#' variable must be a factor variable.
348		#'
349		#' When optimizer is not set, first the default optimizer
350		#' (`L-BFGS-B`) is used to fit the model. If that converges, this is returned.
351		#' If not, the other available optimizers from [h_get_optimizers()],
352		#' including `BFGS`, `CG` and `nlminb` are
353		#' tried (in parallel if `n_cores` is set and not on Windows).
354		#' If none of the optimizers converge, then the function fails. Otherwise
355		#' the best fit is returned.
356		#'
357		#' Note that fine-grained control specifications can either be passed directly
358		#' to the `mmrm` function, or via the `control` argument for bundling together
359		#' with the [mmrm_control()] function. Both cannot be used together, since
360		#' this would delete the arguments passed via `mmrm`.
361		#'
362		#' @return An `mmrm` object.
363		#'
364		#' @note The `mmrm` object is also an `mmrm_fit` and an `mmrm_tmb` object,
365		#' therefore corresponding methods also work (see [`mmrm_tmb_methods`]).
366		#'
367		#' Additional contents depend on the choice of the adjustment `method`:
368		#' - If Satterthwaite adjustment is used, the Jacobian information `jac_list`
369		#' is included.
370		#' - If Kenward-Roger adjustment is used, `kr_comp` contains necessary
371		#' components and `beta_vcov_adj` includes the adjusted coefficients covariance
372		#' matrix.
373		#'
374		#' Use of the package `emmeans` is supported, see [`emmeans_support`].
375		#'
376		#' NA values are always omitted regardless of `na.action` setting.
377		#'
378		#' When the number of visit levels is large, it usually requires large memory to create the
379		#' covariance matrix. By default, the maximum allowed visit levels is 100, and if there are more
380		#' visit levels, a confirmation is needed if run interactively.
381		#' You can use `options(mmrm.max_visits = <target>)` to increase the maximum allowed number of visit
382		#' levels. In non-interactive sessions the confirmation is not raised and will directly give you an error if
383		#' the number of visit levels exceeds the maximum.
384		#'
385		#' @export
386		#'
387		#' @examples
388		#' fit <- mmrm(
389		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
390		#' data = fev_data
391		#' )
392		#'
393		#' # Direct specification of control details:
394		#' fit <- mmrm(
395		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
396		#' data = fev_data,
397		#' weights = fev_data$WEIGHTS,
398		#' method = "Kenward-Roger"
399		#' )
400		#'
401		#' # Alternative specification via control argument (but you cannot mix the
402		#' # two approaches):
403		#' fit <- mmrm(
404		#' formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID),
405		#' data = fev_data,
406		#' control = mmrm_control(method = "Kenward-Roger")
407		#' )
408		mmrm <- function(formula,
409		data,
410		weights = NULL,
411		covariance = NULL,
412		reml = TRUE,
413		control = mmrm_control(...),
414		...) {
415	175x	assert_false(!missing(control) && !missing(...))
416	174x	assert_class(control, "mmrm_control")
417	169x	assert_list(control$optimizers, min.len = 1)
418
419	169x	if (control$method %in% c("Kenward-Roger", "Kenward-Roger-Linear") && !reml) {
420	!	stop("Kenward-Roger only works for REML")
421		}
422	169x	h_valid_formula(formula)
423	168x	covariance <- h_reconcile_cov_struct(formula, covariance)
424	167x	formula_parts <- h_mmrm_tmb_formula_parts(formula, covariance)
425	167x	h_tmb_warn_non_deterministic()
426
427	167x	if (!missing(data)) {
428	166x	attr(data, which = "dataname") <- toString(match.call()$data)
429		} else {
430		# na.action set to na.pass to allow data to be full; will be futher trimmed later
431	1x	data <- model.frame(formula_parts$full_formula, na.action = "na.pass")
432		}
433
434	167x	if (is.null(weights)) {
435	151x	weights <- rep(1, nrow(data))
436		} else {
437	16x	attr(weights, which = "dataname") <- deparse(match.call()$weights)
438		}
439	167x	tmb_data <- h_mmrm_tmb_data(
440	167x	formula_parts, data, weights, reml,
441	167x	singular = if (control$accept_singular) "drop" else "error",
442	167x	drop_visit_levels = control$drop_visit_levels,
443	167x	allow_na_response = FALSE
444		)
445	167x	fit <- structure("", class = "try-error")
446	167x	names_all_optimizers <- names(control$optimizers)
447	167x	while (is(fit, "try-error") && length(control$optimizers) > 0) {
448	171x	fit <- fit_single_optimizer(
449	171x	tmb_data = tmb_data,
450	171x	formula_parts = formula_parts,
451	171x	control = control
452		)
453	168x	if (is(fit, "try-error")) {
454	6x	warning(paste0(
455	6x	"Divergence with optimizer ", names(control$optimizers[1L]), " due to problems: ",
456	6x	toString(attr(fit, "divergence"))
457		))
458		}
459	168x	control$optimizers <- control$optimizers[-1]
460		}
461	164x	if (!attr(fit, "converged")) {
462	7x	more_optimizers <- length(control$optimizers) >= 1L
463	7x	if (more_optimizers) {
464	5x	fit <- refit_multiple_optimizers(
465	5x	fit = fit,
466	5x	control = control
467		)
468		} else {
469	2x	all_problems <- unlist(
470	2x	attributes(fit)[c("errors", "warnings")],
471	2x	use.names = FALSE
472		)
473	2x	stop(paste0(
474	2x	"Chosen optimizers '", toString(names_all_optimizers), "' led to problems during model fit:\n",
475	2x	paste(paste0(seq_along(all_problems), ") ", all_problems), collapse = ";\n"), "\n",
476	2x	"Consider trying multiple or different optimizers."
477		))
478		}
479		}
480	161x	fit_msg <- attr(fit, "messages")
481	161x	if (!is.null(fit_msg)) {
482	!	message(paste(fit_msg, collapse = "\n"))
483		}
484	161x	fit$call <- match.call()
485	161x	fit$call$formula <- formula
486	161x	fit$method <- control$method
487	161x	fit$vcov <- control$vcov
488	161x	if (control$vcov %in% c("Kenward-Roger", "Kenward-Roger-Linear")) {
489	47x	fit$kr_comp <- h_get_kr_comp(fit$tmb_data, fit$theta_est)
490	47x	fit$beta_vcov_adj <- h_var_adj(
491	47x	v = fit$beta_vcov,
492	47x	w = component(fit, "theta_vcov"),
493	47x	p = fit$kr_comp$P,
494	47x	q = fit$kr_comp$Q,
495	47x	r = fit$kr_comp$R,
496	47x	linear = (control$vcov == "Kenward-Roger-Linear")
497		)
498	114x	} else if (control$vcov %in% c("Empirical", "Empirical-Bias-Reduced", "Empirical-Jackknife")) {
499	31x	empirical_comp <- h_get_empirical(
500	31x	fit$tmb_data, fit$theta_est, fit$beta_est, fit$beta_vcov, control$vcov
501		)
502	31x	fit$beta_vcov_adj <- empirical_comp$cov
503	31x	fit$empirical_df_mat <- empirical_comp$df_mat
504	31x	fit$score_per_subject <- empirical_comp$score_per_subject
505	31x	dimnames(fit$beta_vcov_adj) <- dimnames(fit$beta_vcov)
506	83x	} else if (identical(control$vcov, "Asymptotic")) {
507		# Note that we only need the Jacobian list under Asymptotic covariance method,
508		# cf. the Satterthwaite vignette.
509	83x	if (identical(fit$method, "Satterthwaite")) {
510	81x	fit$jac_list <- h_jac_list(fit$tmb_data, fit$theta_est, fit$beta_vcov)
511		}
512		} else {
513	!	stop("Unrecognized coefficent variance-covariance method!")
514		}
515
516	161x	class(fit) <- c("mmrm", class(fit))
517	161x	fit
518		}

1		#' Register `mmrm` For Use With `car::Anova`
2		#'
3		#' @inheritParams base::requireNamespace
4		#' @return A logical value indicating whether registration was successful.
5		#'
6		#' @keywords internal
7		car_add_mmrm <- function(quietly = FALSE) {
8	1x	if (!requireNamespace("car", quietly = quietly)) {
9	!	return(FALSE)
10		}
11	1x	envir <- asNamespace("mmrm")
12	1x	h_register_s3("car", "Anova", "mmrm", envir)
13	1x	TRUE
14		}
15
16
17		#' Obtain Contrast for Specified Effect
18		#'
19		#' This is support function to obtain contrast matrix for type II/III testing.
20		#'
21		#' @param object (`mmrm`)\cr the fitted MMRM.
22		#' @param effect (`string`) the name of the effect.
23		#' @param type (`string`) type of test, "II", "III", '2', or '3'.
24		#' @param tol (`numeric`) threshold blow which values are treated as 0.
25		#'
26		#' @return A `matrix` of the contrast.
27		#'
28		#' @keywords internal
29		h_get_contrast <- function(object, effect, type = c("II", "III", "2", "3"), tol = sqrt(.Machine$double.eps)) {
30	45x	assert_class(object, "mmrm")
31	45x	assert_string(effect)
32	45x	assert_double(tol, finite = TRUE, len = 1L)
33	45x	type <- match.arg(type)
34	45x	mx <- component(object, "x_matrix")
35	45x	asg <- attr(mx, "assign")
36	45x	formula <- object$formula_parts$model_formula
37	45x	tms <- terms(formula)
38	45x	fcts <- attr(tms, "factors")[-1L, , drop = FALSE] # Discard the response.
39	45x	ods <- attr(tms, "order")
40	45x	assert_subset(effect, colnames(fcts))
41	45x	idx <- which(effect == colnames(fcts))
42	45x	cols <- which(asg == idx)
43	45x	xlev <- component(object, "xlev")
44	45x	contains_intercept <- (!0 %in% asg) && h_first_contain_categorical(effect, fcts, names(xlev))
45	45x	coef_rows <- length(cols) - as.integer(contains_intercept)
46	45x	l_mx <- matrix(0, nrow = coef_rows, ncol = length(asg))
47	45x	if (coef_rows == 0L) {
48	1x	return(l_mx)
49		}
50	44x	if (contains_intercept) {
51	4x	l_mx[, cols] <- cbind(-1, diag(rep(1, coef_rows)))
52		} else {
53	40x	l_mx[, cols] <- diag(rep(1, coef_rows))
54		}
55	44x	for (i in setdiff(seq_len(ncol(fcts)), idx)) {
56	120x	additional_vars <- names(which(fcts[, i] > fcts[, idx]))
57	120x	additional_numeric <- any(!additional_vars %in% names(xlev))
58	120x	current_col <- which(asg == i)
59	120x	if (ods[i] >= ods[idx] && all(fcts[, i] >= fcts[, idx]) && !additional_numeric) {
60	24x	sub_mat <- switch(type,
61	24x	"2" = ,
62	24x	"II" = {
63	8x	x1 <- mx[, cols, drop = FALSE]
64	8x	x0 <- mx[, -c(cols, current_col), drop = FALSE]
65	8x	x2 <- mx[, current_col, drop = FALSE]
66	8x	m <- diag(rep(1, nrow(x0))) - x0 %% solve(t(x0) %% x0) %*% t(x0)
67	8x	ret <- solve(t(x1) %% m %% x1) %% t(x1) %% m %*% x2
68	8x	if (contains_intercept) {
69	1x	ret[-1, ] - ret[1, ]
70		} else {
71	7x	ret
72		}
73		},
74	24x	"3" = ,
75	24x	"III" = {
76	16x	lvls <- h_obtain_lvls(effect, additional_vars, xlev)
77	16x	t_levels <- lvls$total
78	16x	nms_base <- colnames(mx)[cols]
79	16x	nms <- colnames(mx)[current_col]
80	16x	nms_base_split <- strsplit(nms_base, ":")
81	16x	nms_split <- strsplit(nms, ":")
82	16x	base_idx <- h_get_index(nms_split, nms_base_split)
83	16x	mt <- l_mx[, cols, drop = FALSE] / t_levels
84	16x	ret <- mt[, base_idx, drop = FALSE]
85		# if there is extra levels, replace it with -1/t_levels
86	16x	ret[is.na(ret)] <- -1 / t_levels
87	16x	ret
88		}
89		)
90	24x	l_mx[, current_col] <- sub_mat
91		}
92		}
93	44x	l_mx[abs(l_mx) < tol] <- 0
94	44x	l_mx
95		}
96
97		#' Conduct type II/III hypothesis testing on the MMRM fit results.
98		#'
99		#' @param mod (`mmrm`)\cr the fitted MMRM.
100		#' @param ... not used.
101		#' @inheritParams h_get_contrast
102		#'
103		#' @details
104		#' `Anova` will return `anova` object with one row per variable and columns
105		#' `Num Df`(numerator degrees of freedom), `Denom Df`(denominator degrees of freedom),
106		#' `F Statistic` and `Pr(>=F)`.
107		#'
108		#' @keywords internal
109		# Please do not load `car` and then create the documentation. The Rd file will be different.
110		Anova.mmrm <- function(mod, type = c("II", "III", "2", "3"), tol = sqrt(.Machine$double.eps), ...) { # nolint
111	9x	assert_double(tol, finite = TRUE, len = 1L)
112	9x	type <- match.arg(type)
113	9x	vars <- colnames(attr(terms(mod$formula_parts$model_formula), "factors"))
114	9x	ret <- lapply(
115	9x	vars,
116	9x	function(x) df_md(mod, h_get_contrast(mod, x, type, tol))
117		)
118	9x	ret_df <- do.call(rbind.data.frame, ret)
119	9x	row.names(ret_df) <- vars
120	9x	colnames(ret_df) <- c("Num Df", "Denom Df", "F Statistic", "Pr(>=F)")
121	9x	class(ret_df) <- c("anova", "data.frame")
122	9x	attr(ret_df, "heading") <- sprintf(
123	9x	"Analysis of Fixed Effect Table (Type %s F tests)",
124	9x	switch(type,
125	9x	"2" = ,
126	9x	"II" = "II",
127	9x	"3" = ,
128	9x	"III" = "III"
129		)
130		)
131	9x	ret_df
132		}
133
134
135		#' Obtain Levels Prior and Posterior
136		#' @param var (`string`) name of the effect.
137		#' @param additional_vars (`character`) names of additional variables.
138		#' @param xlev (`list`) named list of character levels.
139		#' @param factors (`matrix`) the factor matrix.
140		#' @keywords internal
141		h_obtain_lvls <- function(var, additional_vars, xlev, factors) {
142	18x	assert_string(var)
143	18x	assert_character(additional_vars)
144	18x	assert_list(xlev, types = "character")
145	18x	nms <- names(xlev)
146	18x	assert_subset(additional_vars, nms)
147	18x	if (var %in% nms) {
148	14x	prior_vars <- intersect(nms[seq_len(match(var, nms) - 1)], additional_vars)
149	14x	prior_lvls <- vapply(xlev[prior_vars], length, FUN.VALUE = 1L)
150	14x	post_vars <- intersect(nms[seq(match(var, nms) + 1, length(nms))], additional_vars)
151	14x	post_lvls <- vapply(xlev[post_vars], length, FUN.VALUE = 1L)
152	14x	total_lvls <- prod(prior_lvls) * prod(post_lvls)
153		} else {
154	4x	prior_lvls <- vapply(xlev[additional_vars], length, FUN.VALUE = 1L)
155	4x	post_lvls <- 2L
156	4x	total_lvls <- prod(prior_lvls)
157		}
158	18x	list(
159	18x	prior = prior_lvls,
160	18x	post = post_lvls,
161	18x	total = total_lvls
162		)
163		}
164
165		#' Check if the Effect is the First Categorical Effect
166		#' @param effect (`string`) name of the effect.
167		#' @param categorical (`character`) names of the categorical values.
168		#' @param factors (`matrix`) the factor matrix.
169		#' @keywords internal
170		h_first_contain_categorical <- function(effect, factors, categorical) {
171	9x	assert_string(effect)
172	9x	assert_matrix(factors)
173	9x	assert_character(categorical)
174	9x	mt <- match(effect, colnames(factors))
175	9x	varnms <- row.names(factors)
176		# if the effect is not categorical in any value, return FALSE
177	9x	if (!any(varnms[factors[, mt] > 0] %in% categorical)) {
178	2x	return(FALSE)
179		}
180		# keep only categorical rows that is in front of the current factor
181	7x	factors <- factors[row.names(factors) %in% categorical, seq_len(mt - 1L), drop = FALSE]
182		# if previous cols are all numerical, return TRUE
183	7x	if (ncol(factors) < 1L) {
184	4x	return(TRUE)
185		}
186	3x	col_ind <- apply(factors, 2, prod)
187		# if any of the previous cols are categorical, return FALSE
188	3x	return(!any(col_ind > 0))
189		}
190
191		#' Test if the First Vector is Subset of the Second Vector
192		#' @param x (`vector`) the first list.
193		#' @param y (`vector`) the second list.
194		#' @keywords internal
195		h_get_index <- function(x, y) {
196	18x	assert_list(x)
197	18x	assert_list(y)
198	18x	vapply(
199	18x	x,
200	18x	\(i) {
201	68x	r <- vapply(y, \(j) test_subset(j, i), FUN.VALUE = TRUE)
202	68x	if (sum(r) == 1L) {
203	65x	which(r)
204		} else {
205	18x	NA_integer_
206		}
207		},
208	18x	FUN.VALUE = 1L
209		)
210		}

1		#' Obtain List of Jacobian Matrix Entries for Covariance Matrix
2		#'
3		#' @description Obtain the Jacobian matrices given the covariance function and variance parameters.
4		#'
5		#' @param tmb_data (`mmrm_tmb_data`)\cr produced by [h_mmrm_tmb_data()].
6		#' @param theta_est (`numeric`)\cr variance parameters point estimate.
7		#' @param beta_vcov (`matrix`)\cr vairance covariance matrix of coefficients.
8		#'
9		#' @return List with one element per variance parameter containing a matrix
10		#' of the same dimensions as the covariance matrix. The values are the derivatives
11		#' with regards to this variance parameter.
12		#'
13		#' @keywords internal
14		h_jac_list <- function(tmb_data,
15		theta_est,
16		beta_vcov) {
17	82x	assert_class(tmb_data, "mmrm_tmb_data")
18	82x	assert_numeric(theta_est)
19	82x	assert_matrix(beta_vcov)
20	82x	.Call(`_mmrm_get_jacobian`, PACKAGE = "mmrm", tmb_data, theta_est, beta_vcov)
21		}
22
23		#' Quadratic Form Calculations
24		#'
25		#' @description These helpers are mainly for easier readability and slightly better efficiency
26		#' of the quadratic forms used in the Satterthwaite calculations.
27		#'
28		#' @param center (`matrix`)\cr square numeric matrix with the same dimensions as
29		#' `x` as the center of the quadratic form.
30		#'
31		#' @name h_quad_form
32		NULL
33
34		#' @describeIn h_quad_form calculates the number `vec %% center %% t(vec)`
35		#' as a numeric (not a matrix).
36		#'
37		#' @param vec (`numeric`)\cr interpreted as a row vector.
38		#'
39		#' @keywords internal
40		h_quad_form_vec <- function(vec, center) {
41	5607x	vec <- as.vector(vec)
42	5607x	assert_numeric(vec, any.missing = FALSE)
43	5607x	assert_matrix(
44	5607x	center,
45	5607x	mode = "numeric",
46	5607x	any.missing = FALSE,
47	5607x	nrows = length(vec),
48	5607x	ncols = length(vec)
49		)
50
51	5607x	sum(vec * (center %*% vec))
52		}
53
54		#' @describeIn h_quad_form calculates the quadratic form `mat %% center %% t(mat)`
55		#' as a matrix, the result is square and has dimensions identical to the number
56		#' of rows in `mat`.
57		#'
58		#' @param mat (`matrix`)\cr numeric matrix to be multiplied left and right of
59		#' `center`, therefore needs to have as many columns as there are rows and columns
60		#' in `center`.
61		#'
62		#' @keywords internal
63		h_quad_form_mat <- function(mat, center) {
64	119x	assert_matrix(mat, mode = "numeric", any.missing = FALSE, min.cols = 1L)
65	119x	assert_matrix(
66	119x	center,
67	119x	mode = "numeric",
68	119x	any.missing = FALSE,
69	119x	nrows = ncol(center),
70	119x	ncols = ncol(center)
71		)
72	119x	mat %*% tcrossprod(center, mat)
73		}
74
75		#' Computation of a Gradient Given Jacobian and Contrast Vector
76		#'
77		#' @description Computes the gradient of a linear combination of `beta` given the Jacobian matrix and
78		#' variance parameters.
79		#'
80		#' @param jac_list (`list`)\cr Jacobian list produced e.g. by [h_jac_list()].
81		#' @param contrast (`numeric`)\cr contrast vector, which needs to have the
82		#' same number of elements as there are rows and columns in each element of
83		#' `jac_list`.
84		#'
85		#' @return Numeric vector which contains the quadratic forms of each element of
86		#' `jac_list` with the `contrast` vector.
87		#'
88		#' @keywords internal
89		h_gradient <- function(jac_list, contrast) {
90	491x	assert_list(jac_list)
91	491x	assert_numeric(contrast)
92
93	491x	vapply(
94	491x	jac_list,
95	491x	h_quad_form_vec,
96	491x	vec = contrast,
97	491x	numeric(1L)
98		)
99		}
100
101		#' Calculation of Satterthwaite Degrees of Freedom for One-Dimensional Contrast
102		#'
103		#' @description Used in [df_1d()] if method is
104		#' "Satterthwaite".
105		#'
106		#' @param object (`mmrm`)\cr the MMRM fit.
107		#' @param contrast (`numeric`)\cr contrast vector. Note that this should not include
108		#' elements for singular coefficient estimates, i.e. only refer to the
109		#' actually estimated coefficients.
110		#'
111		#' @return List with `est`, `se`, `df`, `t_stat` and `p_val`.
112		#' @keywords internal
113		h_df_1d_sat <- function(object, contrast) {
114	456x	assert_class(object, "mmrm")
115	456x	contrast <- as.numeric(contrast)
116	456x	assert_numeric(contrast, len = length(component(object, "beta_est")))
117
118	456x	df <- if (identical(object$vcov, "Asymptotic")) {
119	444x	grad <- h_gradient(component(object, "jac_list"), contrast)
120	444x	v_num <- 2 * h_quad_form_vec(contrast, component(object, "beta_vcov"))^2
121	444x	v_denom <- h_quad_form_vec(grad, component(object, "theta_vcov"))
122	444x	v_num / v_denom
123	456x	} else if (object$vcov %in% c("Empirical", "Empirical-Jackknife", "Empirical-Bias-Reduced")) {
124	12x	contrast_matrix <- Matrix::.bdiag(rep(list(matrix(contrast, nrow = 1)), component(object, "n_subjects")))
125	12x	contrast_matrix <- as.matrix(contrast_matrix)
126	12x	g_matrix <- h_quad_form_mat(contrast_matrix, object$empirical_df_mat)
127	12x	h_tr(g_matrix)^2 / sum(g_matrix^2)
128		}
129
130	456x	h_test_1d(object, contrast, df)
131		}
132
133		#' Calculating Denominator Degrees of Freedom for the Multi-Dimensional Case
134		#'
135		#' @description Calculates the degrees of freedom for multi-dimensional contrast.
136		#'
137		#' @param t_stat_df (`numeric`)\cr `n` t-statistic derived degrees of freedom.
138		#'
139		#' @return Usually the calculation is returning `2 * E / (E - n)` where
140		#' `E` is the sum of `t / (t - 2)` over all `t_stat_df` values `t`.
141		#'
142		#' @note If the input values are two similar to each other then just the average
143		#' of them is returned. If any of the inputs is not larger than 2 then 2 is
144		#' returned.
145		#'
146		#' @keywords internal
147		h_md_denom_df <- function(t_stat_df) {
148	24x	assert_numeric(t_stat_df, min.len = 1L, lower = .Machine$double.xmin, any.missing = FALSE)
149
150	24x	if (test_scalar(t_stat_df)) {
151	1x	t_stat_df
152	23x	} else if (all(abs(diff(t_stat_df)) < sqrt(.Machine$double.eps))) {
153	1x	mean(t_stat_df)
154	22x	} else if (any(t_stat_df <= 2)) {
155	2x	2
156		} else {
157	20x	e <- sum(t_stat_df / (t_stat_df - 2))
158	20x	2 * e / (e - (length(t_stat_df)))
159		}
160		}
161
162		#' Creating F-Statistic Results from One-Dimensional Contrast
163		#'
164		#' @description Creates multi-dimensional result from one-dimensional contrast from [df_1d()].
165		#'
166		#' @param object (`mmrm`)\cr model fit.
167		#' @param contrast (`numeric`)\cr one-dimensional contrast.
168		#'
169		#' @return The one-dimensional degrees of freedom are calculated and then
170		#' based on that the p-value is calculated.
171		#'
172		#' @keywords internal
173		h_df_md_from_1d <- function(object, contrast) {
174	134x	res_1d <- h_df_1d_sat(object, contrast)
175	134x	list(
176	134x	num_df = 1,
177	134x	denom_df = res_1d$df,
178	134x	f_stat = res_1d$t_stat^2,
179	134x	p_val = stats::pf(q = res_1d$t_stat^2, df1 = 1, df2 = res_1d$df, lower.tail = FALSE)
180		)
181		}
182
183		#' Calculation of Satterthwaite Degrees of Freedom for Multi-Dimensional Contrast
184		#'
185		#' @description Used in [df_md()] if method is "Satterthwaite".
186		#'
187		#' @param object (`mmrm`)\cr the MMRM fit.
188		#' @param contrast (`matrix`)\cr numeric contrast matrix, if given a `numeric`
189		#' then this is coerced to a row vector. Note that this should not include
190		#' elements for singular coefficient estimates, i.e. only refer to the
191		#' actually estimated coefficients.
192		#'
193		#' @return List with `num_df`, `denom_df`, `f_stat` and `p_val` (2-sided p-value).
194		#' @keywords internal
195		h_df_md_sat <- function(object, contrast) {
196	151x	assert_class(object, "mmrm")
197	151x	assert_matrix(contrast, mode = "numeric", any.missing = FALSE, ncols = length(component(object, "beta_est")))
198		# Early return if we are in the one-dimensional case.
199	151x	if (identical(nrow(contrast), 1L)) {
200	132x	return(h_df_md_from_1d(object, contrast))
201		}
202
203	19x	contrast_cov <- h_quad_form_mat(contrast, component(object, "beta_vcov"))
204	19x	eigen_cont_cov <- eigen(contrast_cov)
205	19x	eigen_cont_cov_vctrs <- eigen_cont_cov$vectors
206	19x	eigen_cont_cov_vals <- eigen_cont_cov$values
207
208	19x	eps <- sqrt(.Machine$double.eps)
209	19x	tol <- max(eps * eigen_cont_cov_vals[1], 0)
210	19x	rank_cont_cov <- sum(eigen_cont_cov_vals > tol)
211	19x	assert_number(rank_cont_cov, lower = .Machine$double.xmin)
212	19x	rank_seq <- seq_len(rank_cont_cov)
213	19x	vctrs_cont_prod <- crossprod(eigen_cont_cov_vctrs, contrast)[rank_seq, , drop = FALSE]
214
215		# Early return if rank 1.
216	19x	if (identical(rank_cont_cov, 1L)) {
217	1x	return(h_df_md_from_1d(object, vctrs_cont_prod))
218		}
219
220	18x	t_squared_nums <- drop(vctrs_cont_prod %*% object$beta_est)^2
221	18x	t_squared_denoms <- eigen_cont_cov_vals[rank_seq]
222	18x	t_squared <- t_squared_nums / t_squared_denoms
223	18x	f_stat <- sum(t_squared) / rank_cont_cov
224	18x	t_stat_df_nums <- 2 * eigen_cont_cov_vals^2
225	18x	t_stat_df <- if (identical(object$vcov, "Asymptotic")) {
226	18x	grads_vctrs_cont_prod <- lapply(
227	18x	rank_seq,
228	18x	function(m) h_gradient(component(object, "jac_list"), contrast = vctrs_cont_prod[m, ])
229		)
230	18x	t_stat_df_denoms <- vapply(
231	18x	grads_vctrs_cont_prod,
232	18x	h_quad_form_vec,
233	18x	center = component(object, "theta_vcov"),
234	18x	numeric(1)
235		)
236	18x	t_stat_df_nums / t_stat_df_denoms
237		} else {
238	!	vapply(
239	!	rank_seq,
240	!	function(m) {
241	!	contrast_matrix <- Matrix::.bdiag(
242	!	rep(list(vctrs_cont_prod[m, , drop = FALSE]), component(object, "n_subjects"))
243		)
244	!	contrast_matrix <- as.matrix(contrast_matrix)
245	!	g_matrix <- h_quad_form_mat(contrast_matrix, object$empirical_df_mat)
246	!	h_tr(g_matrix)^2 / sum(g_matrix^2)
247		},
248	!	FUN.VALUE = 0
249		)
250		}
251	18x	denom_df <- h_md_denom_df(t_stat_df)
252
253	18x	list(
254	18x	num_df = rank_cont_cov,
255	18x	denom_df = denom_df,
256	18x	f_stat = f_stat,
257	18x	p_val = stats::pf(q = f_stat, df1 = rank_cont_cov, df2 = denom_df, lower.tail = FALSE)
258		)
259		}

1		#' Methods for `mmrm` Objects
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' @param object (`mmrm`)\cr the fitted MMRM including Jacobian and call etc.
6		#' @param ... not used.
7		#' @return Depends on the method, see Details and Functions.
8		#'
9		#' @details
10		#' While printing the summary of (`mmrm`)\cr object, the following will be displayed:
11		#' 1. Formula. The formula used in the model.
12		#' 2. Data. The data used for analysis, including number of subjects, number of valid observations.
13		#' 3. Covariance. The covariance structure and number of variance parameters.
14		#' 4. Method. Restricted maximum likelihood(REML) or maximum likelihood(ML).
15		#' 5. Model selection criteria. AIC, BIC, log likelihood and deviance.
16		#' 6. Coefficients. Coefficients of the covariates.
17		#' 7. Covariance estimate. The covariance estimate(for each group).
18		#' 1. If the covariance structure is non-spatial, the covariance matrix of all categorical time points available
19		#' in data will be displayed.
20		#' 2. If the covariance structure is spatial, the covariance matrix of two time points with unit distance
21		#' will be displayed.
22		#'
23		#' `confint` is used to obtain the confidence intervals for the coefficients.
24		#' Please note that this is different from the confidence interval of difference
25		#' of least square means from `emmeans`.
26		#'
27		#' @name mmrm_methods
28		#'
29		#' @seealso [`mmrm_tmb_methods`], [`mmrm_tidiers`] for additional methods.
30		#'
31		#' @examples
32		#' formula <- FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT \| USUBJID)
33		#' object <- mmrm(formula, fev_data)
34		NULL
35
36		#' Coefficients Table for MMRM Fit
37		#'
38		#' This is used by [summary.mmrm()] to obtain the coefficients table.
39		#'
40		#' @param object (`mmrm`)\cr model fit.
41		#'
42		#' @return Matrix with one row per coefficient and columns
43		#' `Estimate`, `Std. Error`, `df`, `t value` and `Pr(>\|t\|)`.
44		#'
45		#' @keywords internal
46		h_coef_table <- function(object) {
47	40x	assert_class(object, "mmrm")
48
49	40x	coef_est <- component(object, "beta_est")
50	40x	coef_contrasts <- diag(x = rep(1, length(coef_est)))
51	40x	rownames(coef_contrasts) <- names(coef_est)
52	40x	coef_table <- t(apply(
53	40x	coef_contrasts,
54	40x	MARGIN = 1L,
55	40x	FUN = function(contrast) unlist(df_1d(object, contrast))
56		))
57	40x	assert_names(
58	40x	colnames(coef_table),
59	40x	identical.to = c("est", "se", "df", "t_stat", "p_val")
60		)
61	40x	colnames(coef_table) <- c("Estimate", "Std. Error", "df", "t value", "Pr(>\|t\|)")
62
63	40x	coef_aliased <- component(object, "beta_aliased")
64	40x	if (any(coef_aliased)) {
65	2x	names_coef_na <- names(which(coef_aliased))
66	2x	coef_na_table <- matrix(
67	2x	data = NA,
68	2x	nrow = length(names_coef_na),
69	2x	ncol = ncol(coef_table),
70	2x	dimnames = list(names_coef_na, colnames(coef_table))
71		)
72	2x	coef_table <- rbind(coef_table, coef_na_table)[names(coef_aliased), ]
73		}
74
75	40x	coef_table
76		}
77
78		#' @describeIn mmrm_methods summarizes the MMRM fit results.
79		#' @exportS3Method
80		#' @examples
81		#' # Summary:
82		#' summary(object)
83		summary.mmrm <- function(object, ...) {
84	20x	aic_list <- list(
85	20x	AIC = AIC(object),
86	20x	BIC = BIC(object),
87	20x	logLik = logLik(object),
88	20x	deviance = deviance(object)
89		)
90	20x	coefficients <- h_coef_table(object)
91	20x	call <- stats::getCall(object)
92	20x	components <- component(object, c(
93	20x	"cov_type", "reml", "n_groups", "n_theta",
94	20x	"n_subjects", "n_timepoints", "n_obs",
95	20x	"beta_vcov", "varcor"
96		))
97	20x	components$method <- object$method
98	20x	components$vcov <- object$vcov
99	20x	structure(
100	20x	c(
101	20x	components,
102	20x	list(
103	20x	coefficients = coefficients,
104	20x	n_singular_coefs = sum(component(object, "beta_aliased")),
105	20x	aic_list = aic_list,
106	20x	call = call
107		)
108		),
109	20x	class = "summary.mmrm"
110		)
111		}
112
113		#' Printing MMRM Function Call
114		#'
115		#' This is used in [print.summary.mmrm()].
116		#'
117		#' @param call (`call`)\cr original [mmrm()] function call.
118		#' @param n_obs (`int`)\cr number of observations.
119		#' @param n_subjects (`int`)\cr number of subjects.
120		#' @param n_timepoints (`int`)\cr number of timepoints.
121		#'
122		#' @keywords internal
123		h_print_call <- function(call, n_obs, n_subjects, n_timepoints) {
124	9x	pass <- 0
125	9x	if (!is.null(tmp <- call$formula)) {
126	9x	cat("Formula: ", deparse(tmp), fill = TRUE)
127	9x	rhs <- tmp[[2]]
128	9x	pass <- nchar(deparse(rhs))
129		}
130	9x	if (!is.null(call$data)) {
131	9x	cat(
132	9x	"Data: ", deparse(call$data), "(used", n_obs, "observations from",
133	9x	n_subjects, "subjects with maximum", n_timepoints, "timepoints)",
134	9x	fill = TRUE
135		)
136		}
137		# Display the expression of weights
138	9x	if (!is.null(call$weights)) {
139	4x	cat("Weights: ", deparse(call$weights), fill = TRUE)
140		}
141		}
142
143		#' Printing MMRM Covariance Type
144		#'
145		#' This is used in [print.summary.mmrm()].
146		#'
147		#' @param cov_type (`string`)\cr covariance structure abbreviation.
148		#' @param n_theta (`count`)\cr number of variance parameters.
149		#' @param n_groups (`count`)\cr number of groups.
150		#' @keywords internal
151		h_print_cov <- function(cov_type, n_theta, n_groups) {
152	9x	assert_string(cov_type)
153	9x	assert_count(n_theta, positive = TRUE)
154	9x	assert_count(n_groups, positive = TRUE)
155	9x	cov_definition <- switch(cov_type,
156	9x	us = "unstructured",
157	9x	toep = "Toeplitz",
158	9x	toeph = "heterogeneous Toeplitz",
159	9x	ar1 = "auto-regressive order one",
160	9x	ar1h = "heterogeneous auto-regressive order one",
161	9x	ad = "ante-dependence",
162	9x	adh = "heterogeneous ante-dependence",
163	9x	cs = "compound symmetry",
164	9x	csh = "heterogeneous compound symmetry",
165	9x	sp_exp = "spatial exponential"
166		)
167
168	9x	catstr <- sprintf(
169	9x	"Covariance: %s (%d variance parameters%s)\n",
170	9x	cov_definition,
171	9x	n_theta,
172	9x	ifelse(n_groups == 1L, "", sprintf(" of %d groups", n_groups))
173		)
174	9x	cat(catstr)
175		}
176
177		#' Printing AIC and other Model Fit Criteria
178		#'
179		#' This is used in [print.summary.mmrm()].
180		#'
181		#' @param aic_list (`list`)\cr list as part of from [summary.mmrm()].
182		#' @param digits (`number`)\cr number of decimal places used with [round()].
183		#'
184		#' @keywords internal
185		h_print_aic_list <- function(aic_list,
186		digits = 1) {
187	6x	diag_vals <- round(unlist(aic_list), digits)
188	6x	diag_vals <- format(diag_vals)
189	6x	print(diag_vals, quote = FALSE)
190		}
191
192		#' @describeIn mmrm_methods prints the MMRM fit summary.
193		#' @exportS3Method
194		#' @keywords internal
195		print.summary.mmrm <- function(x,
196		digits = max(3, getOption("digits") - 3),
197		signif.stars = getOption("show.signif.stars"), # nolint
198		...) {
199	5x	cat("mmrm fit\n\n")
200	5x	h_print_call(x$call, x$n_obs, x$n_subjects, x$n_timepoints)
201	5x	h_print_cov(x$cov_type, x$n_theta, x$n_groups)
202	5x	cat("Method: ", x$method, "\n", sep = "")
203	5x	cat("Vcov Method: ", x$vcov, "\n", sep = "")
204	5x	cat("Inference: ")
205	5x	cat(ifelse(x$reml, "REML", "ML"))
206	5x	cat("\n\n")
207	5x	cat("Model selection criteria:\n")
208	5x	h_print_aic_list(x$aic_list)
209	5x	cat("\n")
210	5x	cat("Coefficients: ")
211	5x	if (x$n_singular_coefs > 0) {
212	1x	cat("(", x$n_singular_coefs, " not defined because of singularities)", sep = "")
213		}
214	5x	cat("\n")
215	5x	stats::printCoefmat(
216	5x	x$coefficients,
217	5x	zap.ind = 3,
218	5x	digits = digits,
219	5x	signif.stars = signif.stars
220		)
221	5x	cat("\n")
222	5x	cat("Covariance estimate:\n")
223	5x	if (is.list(x$varcor)) {
224	1x	for (v in names(x$varcor)) {
225	2x	cat(sprintf("Group: %s\n", v))
226	2x	print(round(x$varcor[[v]], digits = digits))
227		}
228		} else {
229	4x	print(round(x$varcor, digits = digits))
230		}
231	5x	cat("\n")
232	5x	invisible(x)
233		}
234
235
236		#' @describeIn mmrm_methods obtain the confidence intervals for the coefficients.
237		#' @exportS3Method
238		#' @examples
239		#' # Confidence Interval:
240		#' confint(object)
241		confint.mmrm <- function(object, parm, level = 0.95, ...) {
242	20x	cf <- coef(object)
243	20x	pnames <- names(cf)
244	20x	if (missing(parm)) {
245	15x	parm <- pnames
246		}
247	20x	assert(
248	20x	check_subset(parm, pnames),
249	20x	check_integerish(parm, lower = 1L, upper = length(cf))
250		)
251	2x	if (is.numeric(parm)) parm <- pnames[parm]
252	18x	assert_number(level, lower = 0, upper = 1)
253	18x	a <- (1 - level) / 2
254	18x	pct <- paste(format(100 * c(a, 1 - a), trim = TRUE, scientific = FALSE, digits = 3), "%")
255	18x	coef_table <- h_coef_table(object)
256	18x	df <- coef_table[parm, "df"]
257	18x	ses <- coef_table[parm, "Std. Error"]
258	18x	fac <- stats::qt(a, df = df)
259	18x	ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, pct))
260	18x	sefac <- ses * fac
261	18x	ci[] <- cf[parm] + c(sefac, -sefac)
262	18x	ci
263		}

1		#' Capture all Output
2		#'
3		#' This function silences all warnings, errors & messages and instead returns a list
4		#' containing the results (if it didn't error), as well as the warnings, errors
5		#' and messages and divergence signals as character vectors.
6		#'
7		#' @param expr (`expression`)\cr to be executed.
8		#' @param remove (`list`)\cr optional list with elements `warnings`, `errors`,
9		#' `messages` which can be character vectors, which will be removed from the
10		#' results if specified.
11		#' @param divergence (`list`)\cr optional list similar as `remove`, but these
12		#' character vectors will be moved to the `divergence` result and signal
13		#' that the fit did not converge.
14		#'
15		#' @return
16		#' A list containing
17		#'
18		#' - `result`: The object returned by `expr` or `list()` if an error was thrown.
19		#' - `warnings`: `NULL` or a character vector if warnings were thrown.
20		#' - `errors`: `NULL` or a string if an error was thrown.
21		#' - `messages`: `NULL` or a character vector if messages were produced.
22		#' - `divergence`: `NULL` or a character vector if divergence messages were caught.
23		#'
24		#' @keywords internal
25		h_record_all_output <- function(expr,
26		remove = list(),
27		divergence = list()) {
28		# Note: We don't need to and cannot assert `expr` here.
29	201x	assert_list(remove, types = "character")
30	201x	assert_list(divergence, types = "character")
31	201x	env <- new.env()
32	201x	result <- withCallingHandlers(
33	201x	withRestarts(
34	201x	expr,
35	201x	muffleStop = function(e) structure(e$message, class = "try-error")
36		),
37	201x	message = function(m) {
38	6x	msg_without_newline <- gsub(m$message, pattern = "\n$", replacement = "")
39	6x	env$message <- c(env$message, msg_without_newline)
40	6x	invokeRestart("muffleMessage")
41		},
42	201x	warning = function(w) {
43	14x	env$warning <- c(env$warning, w$message)
44	14x	invokeRestart("muffleWarning")
45		},
46	201x	error = function(e) {
47	14x	env$error <- c(env$error, e$message)
48	14x	invokeRestart("muffleStop", e)
49		}
50		)
51	201x	list(
52	201x	result = result,
53	201x	warnings = setdiff(env$warning, c(remove$warnings, divergence$warnings)),
54	201x	errors = setdiff(env$error, c(remove$errors, divergence$errors)),
55	201x	messages = setdiff(env$message, c(remove$messages, divergence$messages)),
56	201x	divergence = c(
57	201x	intersect(env$warning, divergence$warnings),
58	201x	intersect(env$error, divergence$errors),
59	201x	intersect(env$message, divergence$messages)
60		)
61		)
62		}
63
64		#' Trace of a Matrix
65		#'
66		#' @description Obtain the trace of a matrix if the matrix is diagonal, otherwise raise an error.
67		#'
68		#' @param x (`matrix`)\cr square matrix input.
69		#'
70		#' @return The trace of the square matrix.
71		#'
72		#' @keywords internal
73		h_tr <- function(x) {
74	1790x	if (nrow(x) != ncol(x)) {
75	1x	stop("x must be square matrix")
76		}
77	1789x	sum(Matrix::diag(x))
78		}
79
80		#' Split Control List
81		#'
82		#' @description Split the [mmrm_control()] object according to its optimizers and use additional arguments
83		#' to replace the elements in the original object.
84		#'
85		#' @param control (`mmrm_control`)\cr object.
86		#' @param ... additional parameters to update the `control` object.
87		#'
88		#' @return A `list` of `mmrm_control` entries.
89		#' @keywords internal
90		h_split_control <- function(control, ...) {
91	8x	assert_class(control, "mmrm_control")
92	8x	l <- length(control$optimizers)
93	8x	lapply(seq_len(l), function(i) {
94	22x	ret <- utils::modifyList(control, list(...))
95	22x	ret$optimizers <- control$optimizers[i]
96	22x	ret
97		})
98		}
99
100		#' Obtain Optimizer according to Optimizer String Value
101		#'
102		#' @description This function creates optimizer functions with arguments.
103		#'
104		#' @param optimizer (`character`)\cr names of built-in optimizers to try, subset
105		#' of "L-BFGS-B", "BFGS", "CG" and "nlminb".
106		#' @param optimizer_fun (`function` or `list` of `function`)\cr alternatively to `optimizer`,
107		#' an optimizer function or a list of optimizer functions can be passed directly here.
108		#' @param optimizer_args (`list`)\cr additional arguments for `optimizer_fun`.
109		#' @param optimizer_control (`list`)\cr passed to argument `control` in `optimizer_fun`.
110		#'
111		#' @details
112		#' If you want to use only the built-in optimizers:
113		#' - `optimizer` is a shortcut to create a list of built-in optimizer functions
114		#' passed to `optimizer_fun`.
115		#' - Allowed are "L-BFGS-B", "BFGS", "CG" (using [stats::optim()] with corresponding method)
116		#' and "nlminb" (using [stats::nlminb()]).
117		#' - Other arguments should go into `optimizer_args`.
118		#'
119		#' If you want to use your own optimizer function:
120		#' - Make sure that there are three arguments: parameter (start value), objective function
121		#' and gradient function are sequentially in the function arguments.
122		#' - If there are other named arguments in front of these, make sure they are correctly
123		#' specified through `optimizer_args`.
124		#' - If the hessian can be used, please make sure its argument name is `hessian` and
125		#' please add attribute `use_hessian = TRUE` to the function,
126		#' using `attr(fun, "use_hessian) <- TRUE`.
127		#'
128		#' @return Named `list` of optimizers created by [h_partial_fun_args()].
129		#'
130		#' @keywords internal
131		h_get_optimizers <- function(optimizer = c("L-BFGS-B", "BFGS", "CG", "nlminb"),
132		optimizer_fun = h_optimizer_fun(optimizer),
133		optimizer_args = list(),
134		optimizer_control = list()) {
135	246x	if ("automatic" %in% optimizer) {
136	1x	lifecycle::deprecate_warn(
137	1x	when = "0.2.0",
138	1x	what = I("\"automatic\" optimizer"),
139	1x	details = "please just omit optimizer argument"
140		)
141	1x	optimizer_fun <- h_optimizer_fun()
142		}
143	246x	assert(
144	246x	test_function(optimizer_fun),
145	246x	test_list(optimizer_fun, types = "function", names = "unique")
146		)
147	246x	if (is.function(optimizer_fun)) {
148	7x	optimizer_fun <- list(custom_optimizer = optimizer_fun)
149		}
150	246x	lapply(optimizer_fun, function(x) {
151	924x	do.call(h_partial_fun_args, c(list(fun = x, control = optimizer_control), optimizer_args))
152		})
153		}
154
155		#' Obtain Optimizer Function with Character
156		#' @description Obtain the optimizer function through the character provided.
157		#' @param optimizer (`character`)\cr vector of optimizers.
158		#'
159		#' @return A (`list`)\cr of optimizer functions generated from [h_partial_fun_args()].
160		#' @keywords internal
161		h_optimizer_fun <- function(optimizer = c("L-BFGS-B", "BFGS", "CG", "nlminb")) {
162	240x	optimizer <- match.arg(optimizer, several.ok = TRUE)
163	240x	lapply(stats::setNames(optimizer, optimizer), function(x) {
164	920x	switch(x,
165	229x	"L-BFGS-B" = h_partial_fun_args(fun = stats::optim, method = x),
166	230x	"BFGS" = h_partial_fun_args(fun = stats::optim, method = x),
167	228x	"CG" = h_partial_fun_args(fun = stats::optim, method = x),
168	233x	"nlminb" = h_partial_fun_args(fun = stats::nlminb, additional_attr = list(use_hessian = TRUE))
169		)
170		})
171		}
172
173		#' Create Partial Functions
174		#' @description Creates partial functions with arguments.
175		#'
176		#' @param fun (`function`)\cr to be wrapped.
177		#' @param ... Additional arguments for `fun`.
178		#' @param additional_attr (`list`)\cr of additional attributes to apply to the result.
179		#'
180		#' @details This function add `args` attribute to the original function,
181		#' and add an extra class `partial` to the function.
182		#' `args` is the argument for the function, and elements in `...` will override the existing
183		#' arguments in attribute `args`. `additional_attr` will override the existing attributes.
184		#'
185		#' @return Object with S3 class `"partial"`, a `function` with `args` attribute (and possibly more
186		#' attributes from `additional_attr`).
187		#' @keywords internal
188		h_partial_fun_args <- function(fun, ..., additional_attr = list()) {
189	1848x	assert_function(fun)
190	1848x	assert_list(additional_attr, names = "unique")
191	1848x	a_args <- list(...)
192	1848x	assert_list(a_args, names = "unique")
193	1848x	args <- attr(fun, "args")
194	1848x	if (is.null(args)) {
195	932x	args <- list()
196		}
197	1848x	do.call(
198	1848x	structure,
199	1848x	args = utils::modifyList(
200	1848x	list(
201	1848x	.Data = fun,
202	1848x	args = utils::modifyList(args, a_args),
203	1848x	class = c("partial", "function")
204		),
205	1848x	additional_attr
206		)
207		)
208		}
209
210		#' Obtain Default Covariance Method
211		#'
212		#' @description Obtain the default covariance method depending on
213		#' the degrees of freedom method used.
214		#'
215		#' @param method (`string`)\cr degrees of freedom method.
216		#'
217		#' @details The default covariance method is different for different degrees of freedom method.
218		#' For "Satterthwaite" or "Between-Within", "Asymptotic" is returned.
219		#' For "Kenward-Roger" only, "Kenward-Roger" is returned.
220		#' For "Residual" only, "Empirical" is returned.
221		#'
222		#' @return String of the default covariance method.
223		#' @keywords internal
224		h_get_cov_default <- function(method = c("Satterthwaite", "Kenward-Roger", "Residual", "Between-Within")) {
225	197x	assert_string(method)
226	197x	method <- match.arg(method)
227	196x	switch(method,
228	1x	"Residual" = "Empirical",
229	158x	"Satterthwaite" = "Asymptotic",
230	35x	"Kenward-Roger" = "Kenward-Roger",
231	2x	"Between-Within" = "Asymptotic"
232		)
233		}
234
235		#' Complete `character` Vector Names From Values
236		#'
237		#' @param x (`character` or `list`)\cr value whose names should be completed
238		#' from element values.
239		#'
240		#' @return A named vector or list.
241		#'
242		#' @keywords internal
243		fill_names <- function(x) {
244	4x	n <- names(x)
245	4x	is_unnamed <- if (is.null(n)) rep_len(TRUE, length(x)) else n == ""
246	4x	names(x)[is_unnamed] <- x[is_unnamed]
247	4x	x
248		}
249
250		#' Drop Items from an Indexible
251		#'
252		#' Drop elements from an indexible object (`vector`, `list`, etc.).
253		#'
254		#' @param x Any object that can be consumed by [seq_along()] and indexed by a
255		#' logical vector of the same length.
256		#' @param n (`integer`)\cr the number of terms to drop.
257		#'
258		#' @return A subset of `x`.
259		#'
260		#' @keywords internal
261		drop_elements <- function(x, n) {
262	819x	x[seq_along(x) > n]
263		}
264
265		#' Ask for Confirmation on Large Visit Levels
266		#'
267		#' @description Ask the user for confirmation if there are too many visit levels
268		#' for non-spatial covariance structure in interactive sessions.
269		#'
270		#' @param x (`numeric`)\cr number of visit levels.
271		#'
272		#' @return Logical value `TRUE`.
273		#' @keywords internal
274		h_confirm_large_levels <- function(x) {
275	297x	assert_count(x)
276	297x	allowed_lvls <- x <= getOption("mmrm.max_visits", 100)
277	297x	if (allowed_lvls) {
278	295x	return(TRUE)
279		}
280	2x	if (!interactive()) {
281	2x	stop("Visit levels too large!", call. = FALSE)
282		}
283	!	proceed <- utils::askYesNo(
284	!	paste(
285	!	"Visit levels is possibly too large.",
286	!	"This requires large memory. Are you sure to continue?",
287	!	collapse = " "
288		)
289		)
290	!	if (!identical(proceed, TRUE)) {
291	!	stop("Visit levels too large!", call. = FALSE)
292		}
293	!	return(TRUE)
294		}
295
296		#' Default Value on NULL
297		#' Return default value when first argument is NULL.
298		#'
299		#' @param x Object.
300		#' @param y Object.
301		#'
302		#' @details If `x` is NULL, returns `y`. Otherwise return `x`.
303		#'
304		#' @keywords internal
305		h_default_value <- function(x, y) {
306	312x	if (is.null(x)) {
307	277x	y
308		} else {
309	35x	x
310		}
311		}
312
313		#' Warn on na.action
314		#' @keywords internal
315		h_warn_na_action <- function() {
316	260x	if (!identical(getOption("na.action"), "na.omit")) {
317	6x	warning("na.action is always set to `na.omit` for `mmrm` fit!")
318		}
319		}
320
321		#' Obtain `na.action` as Function
322		#' @keywords internal
323		h_get_na_action <- function(na_action) {
324	56x	if (is.function(na_action) && identical(methods::formalArgs(na_action), c("object", "..."))) {
325	5x	return(na_action)
326		}
327	51x	if (is.character(na_action) && length(na_action) == 1L) {
328	51x	assert_subset(na_action, c("na.omit", "na.exclude", "na.fail", "na.pass", "na.contiguous"))
329	51x	return(get(na_action, mode = "function", pos = "package:stats"))
330		}
331		}
332
333		#' Validate mmrm Formula
334		#' @param formula (`formula`)\cr to check.
335		#'
336		#' @details In mmrm models, `.` is not allowed as it introduces ambiguity of covariates
337		#' to be used, so it is not allowed to be in formula.
338		#'
339		#' @keywords internal
340		h_valid_formula <- function(formula) {
341	183x	assert_formula(formula)
342	183x	if ("." %in% all.vars(formula)) {
343	2x	stop("`.` is not allowed in mmrm models!")
344		}
345		}
346
347		#' Standard Starting Value
348		#'
349		#' @description Obtain standard start values.
350		#'
351		#' @param cov_type (`string`)\cr name of the covariance structure.
352		#' @param n_visits (`int`)\cr number of visits.
353		#' @param n_groups (`int`)\cr number of groups.
354		#' @param ... not used.
355		#'
356		#' @details
357		#' `std_start` will try to provide variance parameter from identity matrix.
358		#' However, for `ar1` and `ar1h` the corresponding values are not ideal because the
359		#' \eqn{\rho} is usually a positive number thus using 0 as starting value can lead to
360		#' incorrect optimization result, and we use 0.5 as the initial value of \eqn{\rho}.
361		#'
362		#' @return A numeric vector of starting values.
363		#'
364		#' @export
365		std_start <- function(cov_type, n_visits, n_groups, ...) {
366	502x	assert_string(cov_type)
367	502x	assert_subset(cov_type, cov_types(c("abbr", "habbr")))
368	502x	assert_int(n_visits, lower = 1L)
369	502x	assert_int(n_groups, lower = 1L)
370	502x	start_value <- switch(cov_type,
371	502x	us = rep(0, n_visits * (n_visits + 1) / 2),
372	502x	toep = rep(0, n_visits),
373	502x	toeph = rep(0, 2 * n_visits - 1),
374	502x	ar1 = c(0, 0.5),
375	502x	ar1h = c(rep(0, n_visits), 0.5),
376	502x	ad = rep(0, n_visits),
377	502x	adh = rep(0, 2 * n_visits - 1),
378	502x	cs = rep(0, 2),
379	502x	csh = rep(0, n_visits + 1),
380	502x	sp_exp = rep(0, 2)
381		)
382	502x	rep(start_value, n_groups)
383		}
384
385		#' Empirical Starting Value
386		#'
387		#' @description Obtain empirical start value for unstructured covariance
388		#'
389		#' @param data (`data.frame`)\cr data used for model fitting.
390		#' @param model_formula (`formula`)\cr the formula in mmrm model without covariance structure part.
391		#' @param visit_var (`string`)\cr visit variable.
392		#' @param subject_var (`string`)\cr subject id variable.
393		#' @param subject_groups (`factor`)\cr subject group assignment.
394		#' @param ... not used.
395		#'
396		#' @details
397		#' This `emp_start` only works for unstructured covariance structure.
398		#' It uses linear regression to first obtain the coefficients and use the residuals
399		#' to obtain the empirical variance-covariance, and it is then used to obtain the
400		#' starting values.
401		#'
402		#' @note `data` is used instead of `full_frame` because `full_frame` is already
403		#' transformed if model contains transformations, e.g. `log(FEV1) ~ exp(FEV1_BL)` will
404		#' drop `FEV1` and `FEV1_BL` but add `log(FEV1)` and `exp(FEV1_BL)` in `full_frame`.
405		#'
406		#' @return A numeric vector of starting values.
407		#'
408		#' @export
409		em