Package 'mmrm'

Title: Mixed Models for Repeated Measures
Description: Mixed models for repeated measures (MMRM) are a popular choice for analyzing longitudinal continuous outcomes in randomized clinical trials and beyond; see Cnaan, Laird and Slasor (1997) <doi:10.1002/(SICI)1097-0258(19971030)16:20%3C2349::AID-SIM667%3E3.0.CO;2-E> for a tutorial and Mallinckrodt, Lane, Schnell, Peng and Mancuso (2008) <doi:10.1177/009286150804200402> for a review. This package implements MMRM based on the marginal linear model without random effects using Template Model Builder ('TMB') which enables fast and robust model fitting. Users can specify a variety of covariance matrices, weight observations, fit models with restricted or standard maximum likelihood inference, perform hypothesis testing with Satterthwaite or Kenward-Roger adjustment, and extract least square means estimates by using 'emmeans'.
Authors: Daniel Sabanes Bove [aut, cre] , Liming Li [aut], Julia Dedic [aut], Doug Kelkhoff [aut], Kevin Kunzmann [aut], Brian Matthew Lang [aut], Christian Stock [aut], Ya Wang [aut], Craig Gower-Page [ctb], Dan James [aut], Jonathan Sidi [aut], Daniel Leibovitz [aut], Daniel D. Sjoberg [aut] , Lukas A. Widmer [ctb] , Boehringer Ingelheim Ltd. [cph, fnd], Gilead Sciences, Inc. [cph, fnd], F. Hoffmann-La Roche AG [cph, fnd], Merck Sharp & Dohme, Inc. [cph, fnd], AstraZeneca plc [cph, fnd], inferential.biostatistics GmbH [cph, fnd]
Maintainer: Daniel Sabanes Bove <[email protected]>
License: Apache License 2.0
Version: 0.3.14.9001
Built: 2024-11-10 06:21:57 UTC
Source: https://github.com/openpharma/mmrm

Help Index


mmrm Package

Description

mmrm implements mixed models with repeated measures (MMRM) in R.

Author(s)

Maintainer: Daniel Sabanes Bove [email protected] (ORCID)

Authors:

Other contributors:

  • Craig Gower-Page [email protected] [contributor]

  • Lukas A. Widmer (ORCID) [contributor]

  • Boehringer Ingelheim Ltd. [copyright holder, funder]

  • Gilead Sciences, Inc. [copyright holder, funder]

  • F. Hoffmann-La Roche AG [copyright holder, funder]

  • Merck Sharp & Dohme, Inc. [copyright holder, funder]

  • AstraZeneca plc [copyright holder, funder]

  • inferential.biostatistics GmbH [copyright holder, funder]

See Also

Useful links:


Coerce into a Covariance Structure Definition

Description

[Stable]

Usage

as.cov_struct(x, ...)

## S3 method for class 'formula'
as.cov_struct(x, warn_partial = TRUE, ...)

Arguments

x

an object from which to derive a covariance structure. See object specific sections for details.

...

additional arguments unused.

warn_partial

(flag)
whether to emit a warning when parts of the formula are disregarded.

Details

A covariance structure can be parsed from a model definition formula or call. Generally, covariance structures defined using non-standard evaluation take the following form:

type( (visit, )* visit | (group /)? subject )

For example, formulas may include terms such as

us(time | subject)
cp(time | group / subject)
sp_exp(coord1, coord2 | group / subject)

Note that only sp_exp (spatial) covariance structures may provide multiple coordinates, which identify the Euclidean distance between the time points.

Value

A cov_struct() object.

Methods (by class)

  • as.cov_struct(formula): When provided a formula, any specialized functions are assumed to be covariance structure definitions and must follow the form:

    y ~ xs + type( (visit, )* visit | (group /)? subject )
    

    Any component on the right hand side of a formula is considered when searching for a covariance definition.

See Also

Other covariance types: cov_struct(), covariance_types

Examples

# provide a covariance structure as a right-sided formula
as.cov_struct(~ csh(visit | group / subject))

# when part of a full formula, suppress warnings using `warn_partial = FALSE`
as.cov_struct(y ~ x + csh(visit | group / subject), warn_partial = FALSE)

Example Data on BCVA

Description

[Stable]

Usage

bcva_data

Format

A tibble with 10,000 rows and 7 variables:

  • USUBJID: subject ID.

  • VISITN: visit number (numeric).

  • AVISIT: visit number (factor).

  • ARMCD: treatment, TRT or CTL.

  • RACE: 3-category race.

  • BCVA_BL: BCVA at baseline.

  • BCVA_CHG: Change in BCVA at study visits.

Note

Measurements of BCVA (best corrected visual acuity) is a measure of how how many letters a person can read off of an eye chart using corrective lenses or contacts. This a common endpoint in ophthalmology trials.

Source

This is an artificial dataset.


Component Access for mmrm_tmb Objects

Description

[Stable]

Usage

component(
  object,
  name = c("cov_type", "subject_var", "n_theta", "n_subjects", "n_timepoints", "n_obs",
    "beta_vcov", "beta_vcov_complete", "varcor", "formula", "dataset", "n_groups",
    "reml", "convergence", "evaluations", "method", "optimizer", "conv_message", "call",
    "theta_est", "beta_est", "beta_est_complete", "beta_aliased", "x_matrix", "y_vector",
    "neg_log_lik", "jac_list", "theta_vcov", "full_frame", "xlev", "contrasts")
)

Arguments

object

(mmrm_tmb)
the fitted MMRM.

name

(character)
the component(s) to be retrieved.

Details

Available component() names are as follows:

  • call: low-level function call which generated the model.

  • formula: model formula.

  • dataset: data set name.

  • cov_type: covariance structure type.

  • n_theta: number of parameters.

  • n_subjects: number of subjects.

  • n_timepoints: number of modeled time points.

  • n_obs: total number of observations.

  • reml: was REML used (ML was used if FALSE).

  • neg_log_lik: negative log likelihood.

  • convergence: convergence code from optimizer.

  • conv_message: message accompanying the convergence code.

  • evaluations: number of function evaluations for optimization.

  • method: Adjustment method which was used (for mmrm objects), otherwise NULL (for mmrm_tmb objects).

  • beta_vcov: estimated variance-covariance matrix of coefficients (excluding aliased coefficients). When Kenward-Roger/Empirical adjusted coefficients covariance matrix is used, the adjusted covariance matrix is returned (to still obtain the original asymptotic covariance matrix use object$beta_vcov).

  • beta_vcov_complete: estimated variance-covariance matrix including aliased coefficients with entries set to NA.

  • varcor: estimated covariance matrix for residuals. If there are multiple groups, a named list of estimated covariance matrices for residuals will be returned. The names are the group levels.

  • theta_est: estimated variance parameters.

  • beta_est: estimated coefficients (excluding aliased coefficients).

  • beta_est_complete: estimated coefficients including aliased coefficients set to NA.

  • beta_aliased: whether each coefficient was aliased (i.e. cannot be estimated) or not.

  • theta_vcov: estimated variance-covariance matrix of variance parameters.

  • x_matrix: design matrix used (excluding aliased columns).

  • xlev: a named list of character vectors giving the full set of levels to be assumed for each factor.

  • contrasts: a list of contrasts used for each factor.

  • y_vector: response vector used.

  • jac_list: Jacobian, see h_jac_list() for details.

  • full_frame: data.frame with n rows containing all variables needed in the model.

Value

The corresponding component of the object, see details.

See Also

In the lme4 package there is a similar function getME().

Examples

fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)
# Get all available components.
component(fit)
# Get convergence code and message.
component(fit, c("convergence", "conv_message"))
# Get modeled formula as a string.
component(fit, c("formula"))

Define a Covariance Structure

Description

[Stable]

Usage

cov_struct(
  type = cov_types(),
  visits,
  subject,
  group = character(),
  heterogeneous = FALSE
)

Arguments

type

(string)
the name of the covariance structure type to use. For available options, see cov_types(). If a type abbreviation is used that implies heterogeneity (e.g. cph) and no value is provided to heterogeneous, then the heterogeneity is derived from the type name.

visits

(character)
a vector of variable names to use for the longitudinal terms of the covariance structure. Multiple terms are only permitted for the "spatial" covariance type.

subject

(string)
the name of the variable that encodes a subject identifier.

group

(string)
optionally, the name of the variable that encodes a grouping variable for subjects.

heterogeneous

(flag)

Value

A cov_struct object.

See Also

Other covariance types: as.cov_struct(), covariance_types

Examples

cov_struct("csh", "AVISITN", "USUBJID")
cov_struct("spatial", c("VISITA", "VISITB"), group = "GRP", subject = "SBJ")

Covariance Types

Description

[Stable]

Usage

cov_types(
  form = c("name", "abbr", "habbr"),
  filter = c("heterogeneous", "spatial")
)

Arguments

form

(character)
covariance structure type name form. One or more of "name", "abbr" (abbreviation), or "habbr" (heterogeneous abbreviation).

filter

(character)
covariance structure type filter. One or more of "heterogeneous" or "spatial".

Value

A character vector of accepted covariance structure type names and abbreviations.

Abbreviations for Covariance Structures

Common Covariance Structures:

Structure Description Parameters (i,j)(i, j) element
ad Ante-dependence mm σ2k=ij1ρk\sigma^{2}\prod_{k=i}^{j-1}\rho_{k}
adh Heterogeneous ante-dependence 2m12m-1 σiσjk=ij1ρk\sigma_{i}\sigma_{j}\prod_{k=i}^{j-1}\rho_{k}
ar1 First-order auto-regressive 22 σ2ρij\sigma^{2}\rho^{\left \vert {i-j} \right \vert}
ar1h Heterogeneous first-order auto-regressive m+1m+1 σiσjρij\sigma_{i}\sigma_{j}\rho^{\left \vert {i-j} \right \vert}
cs Compound symmetry 22 σ2[ρI(ij)+I(i=j)]\sigma^{2}\left[ \rho I(i \neq j)+I(i=j) \right]
csh Heterogeneous compound symmetry m+1m+1 σiσj[ρI(ij)+I(i=j)]\sigma_{i}\sigma_{j}\left[ \rho I(i \neq j)+I(i=j) \right]
toep Toeplitz mm σij+1\sigma_{\left \vert {i-j} \right \vert +1}
toeph Heterogeneous Toeplitz 2m12m-1 σiσjρij\sigma_{i}\sigma_{j}\rho_{\left \vert {i-j} \right \vert}
us Unstructured m(m+1)/2m(m+1)/2 σij\sigma_{ij}

where ii and jj denote ii-th and jj-th time points, respectively, out of total mm time points, 1i,jm1 \leq i, j \leq m.

Note

The ante-dependence covariance structure in this package refers to homogeneous ante-dependence, while the ante-dependence covariance structure from SAS ⁠PROC MIXED⁠ refers to heterogeneous ante-dependence and the homogeneous version is not available in SAS.

For all non-spatial covariance structures, the time variable must be coded as a factor.

Spatial Covariance structures:

Structure Description Parameters (i,j)(i, j) element
sp_exp spatial exponential 22 σ2ρdij\sigma^{2}\rho^{-d_{ij}}

where dijd_{ij} denotes the Euclidean distance between time points ii and jj.

See Also

Other covariance types: as.cov_struct(), cov_struct()


Calculation of Degrees of Freedom for One-Dimensional Contrast

Description

[Stable] Calculates the estimate, adjusted standard error, degrees of freedom, t statistic and p-value for one-dimensional contrast.

Usage

df_1d(object, contrast)

Arguments

object

(mmrm)
the MMRM fit.

contrast

(numeric)
contrast vector. Note that this should not include elements for singular coefficient estimates, i.e. only refer to the actually estimated coefficients.

Value

List with est, se, df, t_stat and p_val.

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)

Calculation of Degrees of Freedom for Multi-Dimensional Contrast

Description

[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().

Usage

df_md(object, contrast)

Arguments

object

(mmrm)
the MMRM fit.

contrast

(matrix)
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.

Value

List with num_df, denom_df, f_stat and p_val (2-sided p-value).

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)

Support for emmeans

Description

[Stable]

This package includes methods that allow mmrm objects to be used with the emmeans package. emmeans computes estimated marginal means (also called least-square means) for the coefficients of the MMRM. We can also e.g. obtain differences between groups by applying pairs() on the object returned by emmeans::emmeans().

Examples

fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)
if (require(emmeans)) {
  emmeans(fit, ~ ARMCD | AVISIT)
  pairs(emmeans(fit, ~ ARMCD | AVISIT), reverse = TRUE)
}

Empirical Starting Value

Description

Obtain empirical start value for unstructured covariance

Usage

emp_start(data, model_formula, visit_var, subject_var, subject_groups, ...)

Arguments

data

(data.frame)
data used for model fitting.

model_formula

(formula)
the formula in mmrm model without covariance structure part.

visit_var

(string)
visit variable.

subject_var

(string)
subject id variable.

subject_groups

(factor)
subject group assignment.

...

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.

Value

A numeric vector of 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.


Example Data on FEV1

Description

[Stable]

Usage

fev_data

Format

A tibble with 800 rows and 7 variables:

  • USUBJID: subject ID.

  • AVISIT: visit number.

  • ARMCD: treatment, TRT or PBO.

  • RACE: 3-category race.

  • SEX: sex.

  • FEV1_BL: FEV1 at baseline (%).

  • FEV1: FEV1 at study visits.

  • WEIGHT: weighting variable.

  • VISITN: integer order of the visit.

  • VISITN2: coordinates of the visit for distance calculation.

Note

Measurements of FEV1 (forced expired volume in one second) is a measure of how quickly the lungs can be emptied. Low levels of FEV1 may indicate chronic obstructive pulmonary disease (COPD).

Source

This is an artificial dataset.


Low-Level Fitting Function for MMRM

Description

[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().

Usage

fit_mmrm(
  formula,
  data,
  weights,
  reml = TRUE,
  covariance = NULL,
  tmb_data,
  formula_parts,
  control = mmrm_control()
)

Arguments

formula

(formula)
model formula with exactly one special term specifying the visits within subjects, see details.

data

(data.frame)
input data containing the variables used in formula.

weights

(vector)
input vector containing the weights.

reml

(flag)
whether restricted maximum likelihood (REML) estimation is used, otherwise maximum likelihood (ML) is used.

covariance

(cov_struct)
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.

tmb_data

(mmrm_tmb_data)
object.

formula_parts

(mmrm_tmb_formula_parts)
list with formula parts from h_mmrm_tmb_formula_parts().

control

(mmrm_control)
list of control options produced by mmrm_control().

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.

Value

List of class mmrm_tmb, see h_mmrm_tmb_fit() for details. In addition, it contains elements call and optimizer.

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))))

Fitting an MMRM with Single Optimizer

Description

[Stable]

This function helps to fit an MMRM using TMB with a single optimizer, while capturing messages and warnings.

Usage

fit_single_optimizer(
  formula,
  data,
  weights,
  reml = TRUE,
  covariance = NULL,
  tmb_data,
  formula_parts,
  ...,
  control = mmrm_control(...)
)

Arguments

formula

(formula)
the model formula, see details.

data

(data)
the data to be used for the model.

weights

(vector)
an optional vector of weights to be used in the fitting process. Should be NULL or a numeric vector.

reml

(flag)
whether restricted maximum likelihood (REML) estimation is used, otherwise maximum likelihood (ML) is used.

covariance

(cov_struct)
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.

tmb_data

(mmrm_tmb_data)
object.

formula_parts

(mmrm_tmb_formula_parts)
object.

...

Additional arguments to pass to mmrm_control().

control

(mmrm_control)
object.

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.

Value

The mmrm_fit object, with additional attributes containing warnings, messages, optimizer used and convergence status in addition to the mmrm_tmb contents.

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")

Format Covariance Structure Object

Description

Format Covariance Structure Object

Usage

## S3 method for class 'cov_struct'
format(x, ...)

Arguments

x

(cov_struct)
a covariance structure object.

...

Additional arguments unused.

Value

A formatted string for x.


Fit an MMRM

Description

[Stable]

This is the main function fitting the MMRM.

Usage

mmrm(
  formula,
  data,
  weights = NULL,
  covariance = NULL,
  reml = TRUE,
  control = mmrm_control(...),
  ...
)

Arguments

formula

(formula)
the model formula, see details.

data

(data)
the data to be used for the model.

weights

(vector)
an optional vector of weights to be used in the fitting process. Should be NULL or a numeric vector.

covariance

(cov_struct)
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.

reml

(flag)
whether restricted maximum likelihood (REML) estimation is used, otherwise maximum likelihood (ML) is used.

control

(mmrm_control)
fine-grained fitting specifications list created with mmrm_control().

...

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.

Value

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.

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")
)

Control Parameters for Fitting an MMRM

Description

[Stable] Fine-grained specification of the MMRM fit details is possible using this control function.

Usage

mmrm_control(
  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(...)
)

Arguments

n_cores

(count)
number of cores to be used.

method

(string)
adjustment method for degrees of freedom.

vcov

(string)
coefficients covariance matrix adjustment method.

start

(NULL, numeric or function)
optional start values for variance parameters. See details for more information.

accept_singular

(flag)
whether singular design matrices are reduced to full rank automatically and additional coefficient estimates will be missing.

drop_visit_levels

(flag)
whether to drop levels for visit variable, if visit variable is a factor, see details.

...

additional arguments passed to h_get_optimizers().

optimizers

(list)
optimizer specification, created with h_get_optimizers().

Details

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.

Value

List of class mmrm_control with the control parameters.

Examples

mmrm_control(
  optimizer_fun = stats::optim,
  optimizer_args = list(method = "L-BFGS-B")
)

Tidying Methods for mmrm Objects

Description

[Stable]

These methods tidy the estimates from an mmrm object into a summary.

Usage

## S3 method for class 'mmrm'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'mmrm'
glance(x, ...)

## S3 method for class 'mmrm'
augment(
  x,
  newdata = NULL,
  interval = c("none", "confidence", "prediction"),
  se_fit = (interval != "none"),
  type.residuals = c("response", "pearson", "normalized"),
  ...
)

Arguments

x

(mmrm)
fitted model.

conf.int

(flag)
if TRUE columns for the lower (conf.low) and upper bounds (conf.high) of coefficient estimates are included.

conf.level

(number)
defines the range of the optional confidence internal.

...

only used by augment() to pass arguments to the predict.mmrm_tmb() method.

newdata

(data.frame or NULL)
optional new data frame.

interval

(string)
type of interval calculation.

se_fit

(flag)
whether to return standard errors of fit.

type.residuals

(string)
passed on to residuals.mmrm_tmb().

Functions

  • tidy(mmrm): derives tidy tibble from an mmrm object.

  • glance(mmrm): derives glance tibble from an mmrm object.

  • augment(mmrm): derives augment tibble from an mmrm object.

See Also

mmrm_methods, mmrm_tmb_methods for additional methods.

Examples

fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)
# Applying tidy method to return summary table of covariate estimates.
fit |> tidy()
fit |> tidy(conf.int = TRUE, conf.level = 0.9)
# Applying glance method to return summary table of goodness of fit statistics.
fit |> glance()
# Applying augment method to return merged `tibble` of model data, fitted and residuals.
fit |> augment()
fit |> augment(interval = "confidence")
fit |> augment(type.residuals = "pearson")

Methods for mmrm_tmb Objects

Description

[Stable]

Usage

## S3 method for class 'mmrm_tmb'
coef(object, complete = TRUE, ...)

## S3 method for class 'mmrm_tmb'
fitted(object, ...)

## S3 method for class 'mmrm_tmb'
predict(
  object,
  newdata,
  se.fit = FALSE,
  interval = c("none", "confidence", "prediction"),
  level = 0.95,
  nsim = 1000L,
  conditional = FALSE,
  ...
)

## S3 method for class 'mmrm_tmb'
model.frame(
  formula,
  data,
  include = c("subject_var", "visit_var", "group_var", "response_var"),
  full,
  na.action = "na.omit",
  ...
)

## S3 method for class 'mmrm_tmb'
model.matrix(object, data, use_response = TRUE, ...)

## S3 method for class 'mmrm_tmb'
terms(x, include = "response_var", ...)

## S3 method for class 'mmrm_tmb'
logLik(object, ...)

## S3 method for class 'mmrm_tmb'
formula(x, ...)

## S3 method for class 'mmrm_tmb'
vcov(object, complete = TRUE, ...)

## S3 method for class 'mmrm_tmb'
VarCorr(x, sigma = NA, ...)

## S3 method for class 'mmrm_tmb'
deviance(object, ...)

## S3 method for class 'mmrm_tmb'
AIC(object, corrected = FALSE, ..., k = 2)

## S3 method for class 'mmrm_tmb'
BIC(object, ...)

## S3 method for class 'mmrm_tmb'
print(x, ...)

## S3 method for class 'mmrm_tmb'
residuals(object, type = c("response", "pearson", "normalized"), ...)

## S3 method for class 'mmrm_tmb'
simulate(
  object,
  nsim = 1,
  seed = NULL,
  newdata,
  ...,
  method = c("conditional", "marginal")
)

Arguments

object

(mmrm_tmb)
the fitted MMRM object.

complete

(flag)
whether to include potential non-estimable coefficients.

...

mostly not used; Exception is model.matrix() passing ... to the default method.

newdata

(data.frame)
optional new data, otherwise data from object is used.

se.fit

(flag)
indicator if standard errors are required.

interval

(string)
type of interval calculation. Can be abbreviated.

level

(number)
tolerance/confidence level.

nsim

(count)
number of simulations to use.

conditional

(flag)
indicator if the prediction is conditional on the observation or not.

formula

(mmrm_tmb)
same as object.

data

(data.frame)
object in which to construct the frame.

include

(character)
names of variable types to include. Must be NULL or one or more of c("subject_var", "visit_var", "group_var", "response_var").

full

(flag)
indicator whether to return full model frame (deprecated).

na.action

(string)
na action.

use_response

(flag)
whether to use the response for complete rows.

x

(mmrm_tmb)
same as object.

sigma

cannot be used (this parameter does not exist in MMRM).

corrected

(flag)
whether corrected AIC should be calculated.

k

(number)
the penalty per parameter to be used; default k = 2 is the classical AIC.

type

(string)
unscaled (response), pearson or normalized. Default is response, and this is the only type available for use with models with a spatial covariance structure.

seed

unused argument from stats::simulate().

method

(string)
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.

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.

Value

Depends on the method, see Functions.

Functions

  • coef(mmrm_tmb): obtains the estimated coefficients.

  • fitted(mmrm_tmb): obtains the fitted values.

  • predict(mmrm_tmb): 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.

  • model.frame(mmrm_tmb): obtains the model frame.

  • model.matrix(mmrm_tmb): obtains the model matrix.

  • terms(mmrm_tmb): obtains the terms object.

  • logLik(mmrm_tmb): obtains the attained log likelihood value.

  • formula(mmrm_tmb): obtains the used formula.

  • vcov(mmrm_tmb): obtains the variance-covariance matrix estimate for the coefficients.

  • VarCorr(mmrm_tmb): obtains the variance-covariance matrix estimate for the residuals.

  • deviance(mmrm_tmb): 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.

  • AIC(mmrm_tmb): 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 (Hurvich and Tsai 1989; Burnham and Anderson 1998).

  • BIC(mmrm_tmb): obtains the Bayesian Information Criterion, which is using the natural logarithm of the number of subjects for the penalty parameter k.

  • print(mmrm_tmb): prints the object.

  • residuals(mmrm_tmb): to obtain residuals - either unscaled ('response'), 'pearson' or 'normalized'.

  • simulate(mmrm_tmb): 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.

References

  • Hurvich CM, Tsai C (1989). “Regression and time series model selection in small samples.” Biometrika, 76(2), 297–307. doi:10.2307/2336663.

  • Burnham KP, Anderson DR (1998). “Practical use of the information-theoretic approach.” In Model selection and inference, 75–117. Springer. doi:10.1007/978-1-4757-2917-7_3.

  • Gałecki A, Burzykowski T (2013). “Linear mixed-effects model.” In Linear mixed-effects models using R, 245–273. Springer.

See Also

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)))
# Estimated coefficients:
coef(object)
# Fitted values:
fitted(object)
predict(object, newdata = fev_data)
# Model frame:
model.frame(object)
model.frame(object, include = "subject_var")
# Model matrix:
model.matrix(object)
# terms:
terms(object)
terms(object, include = "subject_var")
# Log likelihood given the estimated parameters:
logLik(object)
# Formula which was used:
formula(object)
# Variance-covariance matrix estimate for coefficients:
vcov(object)
# Variance-covariance matrix estimate for residuals:
VarCorr(object)
# REML criterion (twice the negative log likelihood):
deviance(object)
# AIC:
AIC(object)
AIC(object, corrected = TRUE)
# BIC:
BIC(object)
# residuals:
residuals(object, type = "response")
residuals(object, type = "pearson")
residuals(object, type = "normalized")

Print a Covariance Structure Object

Description

Print a Covariance Structure Object

Usage

## S3 method for class 'cov_struct'
print(x, ...)

Arguments

x

(cov_struct)
a covariance structure object.

...

Additional arguments unused.

Value

x invisibly.


Refitting MMRM with Multiple Optimizers

Description

[Stable]

Usage

refit_multiple_optimizers(fit, ..., control = mmrm_control(...))

Arguments

fit

(mmrm_fit)
original model fit from fit_single_optimizer().

...

Additional arguments passed to mmrm_control().

control

(mmrm_control)
object.

Value

The best (in terms of log likelihood) fit which converged.

Note

For Windows, no parallel computations are currently implemented.

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)

Standard Starting Value

Description

Obtain standard start values.

Usage

std_start(cov_type, n_visits, n_groups, ...)

Arguments

cov_type

(string)
name of the covariance structure.

n_visits

(int)
number of visits.

n_groups

(int)
number of groups.

...

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 ρ\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 ρ\rho.

Value

A numeric vector of starting values.