Introduction to `mmrm` • mmrm

Common usage

A minimal call of mmrm(), consisting of only formula and data arguments will produce an object of class mmrm, mmrm_fit, and mmrm_tmb. Here we fit a mmrm model with us (unstructured) covariance structure specified, as well as the defaults of reml = TRUE and control = mmrm_control().

fit <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)

Printing the object will show you output which should be familiar to anyone who has used any popular modeling functions such as stats::lm(), stats::glm(), glmmTMB::glmmTMB(), and lme4::nlmer(). From this print out we see the function call, the data used, the covariance structure with number of variance parameters, as well as the likelihood method, and model deviance achieved. Additionally the user is provided a printout of the estimated coefficients and the model convergence information.

print(fit)
#> mmrm fit
#> 
#> Formula:     FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID)
#> Data:        fev_data (used 537 observations from 197 subjects with maximum 4 
#> timepoints)
#> Covariance:  unstructured (10 variance parameters)
#> Method:      REML
#> Deviance:    3387.373
#> 
#> Coefficients: 
#>                   (Intercept) RACEBlack or African American 
#>                   30.96769899                    1.50464863 
#>                     RACEWhite                      ARMCDTRT 
#>                    5.61309565                    3.77555734 
#>                    AVISITVIS2                    AVISITVIS3 
#>                    4.82858803                   10.33317002 
#>                    AVISITVIS4           ARMCDTRT:AVISITVIS2 
#>                   15.05255715                   -0.01737409 
#>           ARMCDTRT:AVISITVIS3           ARMCDTRT:AVISITVIS4 
#>                   -0.66753189                    0.63094392 
#> 
#> Model Inference Optimization:
#> Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmch

Common customizations

From the high-level mmrm() interface, common changes to the default function call can be specified.

Control Function

For fine control, mmrm_control() is provided. This function allows the user to choose the adjustment method for the degrees of freedom and the coefficients covariance matrix, specify optimization routines, number of cores to be used on Unix systems for trying several optimizers in parallel, provide a vector of starting parameter values, decide the action to be taken when the defined design matrix is singular, not drop unobserved visit levels. For example:

mmrm_control(
  method = "Kenward-Roger",
  optimizer = c("L-BFGS-B", "BFGS"),
  n_cores = 2,
  start = c(0, 1, 1, 0, 1, 0),
  accept_singular = FALSE,
  drop_visit_levels = FALSE
)

Note that this control list can either be passed via the control argument to mmrm, or selected controls can be directly specified in the mmrm call. We will see this below.

REML or ML

Users can specify if REML should be used (default) or if ML should be used in optimization.

fit_ml <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data,
  reml = FALSE
)
fit_ml
#> mmrm fit
#> 
#> Formula:     FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID)
#> Data:        fev_data (used 537 observations from 197 subjects with maximum 4 
#> timepoints)
#> Covariance:  unstructured (10 variance parameters)
#> Method:      ML
#> Deviance:    3397.934
#> 
#> Coefficients: 
#>                   (Intercept) RACEBlack or African American 
#>                    30.9663423                     1.5086851 
#>                     RACEWhite                      ARMCDTRT 
#>                     5.6133151                     3.7761037 
#>                    AVISITVIS2                    AVISITVIS3 
#>                     4.8270155                    10.3353319 
#>                    AVISITVIS4           ARMCDTRT:AVISITVIS2 
#>                    15.0487715                    -0.0156154 
#>           ARMCDTRT:AVISITVIS3           ARMCDTRT:AVISITVIS4 
#>                    -0.6663598                     0.6317222 
#> 
#> Model Inference Optimization:
#> Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmch

Optimizer

Users can specify which optimizer should be used, changing from the default of four optimizers, which starts with L-BFGS-B and proceeds through the other choices if optimization fails to converge. Other choices are BFGS, CG, and nlminb.

L-BFGS-B, BFGS and CG are all implemented with stats::optim() and the Hessian is not used, while nlminb is using stats::nlminb() which in turn uses both the gradient and the Hessian (by default but can be switch off) for the optimization.

fit_opt <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data,
  optimizer = "BFGS"
)
fit_opt
#> mmrm fit
#> 
#> Formula:     FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID)
#> Data:        fev_data (used 537 observations from 197 subjects with maximum 4 
#> timepoints)
#> Covariance:  unstructured (10 variance parameters)
#> Method:      REML
#> Deviance:    3387.373
#> 
#> Coefficients: 
#>                   (Intercept) RACEBlack or African American 
#>                   30.96769004                    1.50467449 
#>                     RACEWhite                      ARMCDTRT 
#>                    5.61310489                    3.77554263 
#>                    AVISITVIS2                    AVISITVIS3 
#>                    4.82858562                   10.33317692 
#>                    AVISITVIS4           ARMCDTRT:AVISITVIS2 
#>                   15.05257072                   -0.01735093 
#>           ARMCDTRT:AVISITVIS3           ARMCDTRT:AVISITVIS4 
#>                   -0.66751927                    0.63095827 
#> 
#> Model Inference Optimization:
#> Converged with code 0 and message:

Covariance Structure

Covariance structures supported by the mmrm are being continuously developed. For a complete list and description please visit the covariance vignette. Below we see the function call for homogeneous compound symmetry (cs).

fit_cs <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + cs(AVISIT | USUBJID),
  data = fev_data,
  reml = FALSE
)
fit_cs
#> mmrm fit
#> 
#> Formula:     FEV1 ~ RACE + ARMCD * AVISIT + cs(AVISIT | USUBJID)
#> Data:        fev_data (used 537 observations from 197 subjects with maximum 4 
#> timepoints)
#> Covariance:  compound symmetry (2 variance parameters)
#> Method:      ML
#> Deviance:    3536.989
#> 
#> Coefficients: 
#>                   (Intercept) RACEBlack or African American 
#>                    31.4207077                     0.5357237 
#>                     RACEWhite                      ARMCDTRT 
#>                     5.4546329                     3.4305212 
#>                    AVISITVIS2                    AVISITVIS3 
#>                     4.8326353                    10.2395076 
#>                    AVISITVIS4           ARMCDTRT:AVISITVIS2 
#>                    15.0672680                     0.2801641 
#>           ARMCDTRT:AVISITVIS3           ARMCDTRT:AVISITVIS4 
#>                    -0.5894964                     0.7939750 
#> 
#> Model Inference Optimization:
#> Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmch

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.

Weighting

Users can perform weighted MMRM by specifying a numeric vector weights with positive values.

fit_wt <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data,
  weights = fev_data$WEIGHT
)
fit_wt
#> mmrm fit
#> 
#> Formula:     FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID)
#> Data:        fev_data (used 537 observations from 197 subjects with maximum 4 
#> timepoints)
#> Weights:     fev_data$WEIGHT
#> Covariance:  unstructured (10 variance parameters)
#> Method:      REML
#> Deviance:    3476.526
#> 
#> Coefficients: 
#>                   (Intercept) RACEBlack or African American 
#>                   31.20065229                    1.18452837 
#>                     RACEWhite                      ARMCDTRT 
#>                    5.36525917                    3.39695951 
#>                    AVISITVIS2                    AVISITVIS3 
#>                    4.85890820                   10.03942420 
#>                    AVISITVIS4           ARMCDTRT:AVISITVIS2 
#>                   14.79354054                    0.03418184 
#>           ARMCDTRT:AVISITVIS3           ARMCDTRT:AVISITVIS4 
#>                    0.01308088                    0.86701567 
#> 
#> Model Inference Optimization:
#> Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmch

Grouped Covariance Structure

Grouped covariance structures are supported by themmrm package. Covariance matrices for each group are identically structured (unstructured, compound symmetry, etc) but the estimates are allowed to vary across groups. We use the form cs(time | group / subject) to specify the group variable.

Here is an example of how we use ARMCD as group variable.

fit_cs <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + cs(AVISIT | ARMCD / USUBJID),
  data = fev_data,
  reml = FALSE
)
VarCorr(fit_cs)
#> $PBO
#>           VIS1      VIS2      VIS3      VIS4
#> VIS1 37.823638  3.601296  3.601296  3.601296
#> VIS2  3.601296 37.823638  3.601296  3.601296
#> VIS3  3.601296  3.601296 37.823638  3.601296
#> VIS4  3.601296  3.601296  3.601296 37.823638
#> 
#> $TRT
#>          VIS1     VIS2     VIS3     VIS4
#> VIS1 49.58110 10.98112 10.98112 10.98112
#> VIS2 10.98112 49.58110 10.98112 10.98112
#> VIS3 10.98112 10.98112 49.58110 10.98112
#> VIS4 10.98112 10.98112 10.98112 49.58110

We can see that the estimated covariance matrices are different in different ARMCD groups.

Adjustment Method

Both Satterthwaite and Kenward-Roger adjustment methods are available. The default is Satterthwaite adjustment of the degrees of freedom. To use Kenward-Roger adjustment of the degrees of freedom as well as the coefficients covariance matrix, use the method argument:

fit_kr <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data,
  method = "Kenward-Roger"
)

Note that this requires reml = TRUE, i.e. Kenward-Roger adjustment is not possible when using maximum likelihood inference. While this adjustment choice is not visible in the print() result of the fitted model (because the initial model fit is not affected by the choice of the adjustment method), looking at the summary we see the method and the correspondingly adjusted standard errors and degrees of freedom:

summary(fit_kr)
#> mmrm fit
#> 
#> Formula:     FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID)
#> Data:        fev_data (used 537 observations from 197 subjects with maximum 4 
#> timepoints)
#> Covariance:  unstructured (10 variance parameters)
#> Method:      Kenward-Roger
#> Inference:   REML
#> 
#> Model selection criteria:
#>      AIC      BIC   logLik deviance 
#>   3407.4   3440.2  -1693.7   3387.4 
#> 
#> Coefficients: 
#>                                Estimate Std. Error        df t value Pr(>|t|)
#> (Intercept)                    30.96770    0.83335 187.91000  37.160  < 2e-16
#> RACEBlack or African American   1.50465    0.62901 169.95000   2.392  0.01784
#> RACEWhite                       5.61310    0.67139 158.87000   8.360 2.98e-14
#> ARMCDTRT                        3.77556    1.07910 146.27000   3.499  0.00062
#> AVISITVIS2                      4.82859    0.80408 143.66000   6.005 1.49e-08
#> AVISITVIS3                     10.33317    0.82303 155.66000  12.555  < 2e-16
#> AVISITVIS4                     15.05256    1.30180 138.39000  11.563  < 2e-16
#> ARMCDTRT:AVISITVIS2            -0.01737    1.13154 138.39000   0.015  0.98777
#> ARMCDTRT:AVISITVIS3            -0.66753    1.18714 158.21000   0.562  0.57470
#> ARMCDTRT:AVISITVIS4             0.63094    1.83319 129.64000   0.344  0.73127
#>                                  
#> (Intercept)                   ***
#> RACEBlack or African American *  
#> RACEWhite                     ***
#> ARMCDTRT                      ***
#> AVISITVIS2                    ***
#> AVISITVIS3                    ***
#> AVISITVIS4                    ***
#> ARMCDTRT:AVISITVIS2              
#> ARMCDTRT:AVISITVIS3              
#> ARMCDTRT:AVISITVIS4              
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Covariance estimate:
#>         VIS1    VIS2    VIS3    VIS4
#> VIS1 40.7335 14.2740  5.1411 13.5288
#> VIS2 14.2740 26.2243  2.6391  7.3219
#> VIS3  5.1411  2.6391 14.9497  1.0341
#> VIS4 13.5288  7.3219  1.0341 95.6006

For one-dimensional contrasts as in the coefficients table above, the degrees of freedom are the same for Kenward-Roger and Satterthwaite. However, Kenward-Roger uses adjusted standard errors, hence the p-values are different.

Note that if you would like to match SAS results for an unstructured covariance model, you can use the linear Kenward-Roger approximation:

fit_kr_lin <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data,
  method = "Kenward-Roger-Linear"
)

This is due to the different parametrization of the unstructured covariance matrix, see the Kenward-Roger vignette for details.

Keeping Unobserved Visits

Sometimes not all possible time points are observed in a given data set. When using a structured covariance matrix, e.g. with auto-regressive structure, then it can be relevant to keep the correct distance between the observed time points.

Consider the following example where we have deliberately removed the VIS3 observations from our initial example data set fev_data to obtain sparse_data. We first fit the model where we do not drop the visit level explicitly, using the drop_visit_levels = FALSE choice. Second we fit the model by default without this option.

sparse_data <- fev_data[fev_data$AVISIT != "VIS3", ]
sparse_result <- mmrm(
  FEV1 ~ RACE + ar1(AVISIT | USUBJID),
  data = sparse_data,
  drop_visit_levels = FALSE
)

dropped_result <- mmrm(
  FEV1 ~ RACE + ar1(AVISIT | USUBJID),
  data = sparse_data
)
#> In AVISIT there are dropped visits: VIS3

We see that we get a message about the dropped visit level by default. Now we can compare the estimated correlation matrices:

cov2cor(VarCorr(sparse_result))
#>            VIS1      VIS2      VIS3       VIS4
#> VIS1 1.00000000 0.4051834 0.1641736 0.06652042
#> VIS2 0.40518341 1.0000000 0.4051834 0.16417360
#> VIS3 0.16417360 0.4051834 1.0000000 0.40518341
#> VIS4 0.06652042 0.1641736 0.4051834 1.00000000
cov2cor(VarCorr(dropped_result))
#>            VIS1      VIS2       VIS4
#> VIS1 1.00000000 0.1468464 0.02156386
#> VIS2 0.14684640 1.0000000 0.14684640
#> VIS4 0.02156386 0.1468464 1.00000000

We see that when using the default, second result, we just drop the VIS3 from the covariance matrix. As a consequence, we model the correlation between VIS2 and VIS4 the same as the correlation between VIS1 and VIS2. Hence we get a smaller correlation estimate here compared to the first result, which includes VIS3 explicitly.

Extraction of model features

Similar to model objects created in other packages, components of mmrm and mmrm_tmb objects can be accessed with standard methods. Additionally, component() is provided to allow deeper and more precise access for those interested in digging through model output. Complete documentation of standard model output methods supported for mmrm_tmb objects can be found at the package website.

Summary

The summary method for mmrm objects provides easy access to frequently needed model components.

fit <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)
fit_summary <- summary(fit)

From this summary object, you can easily retrieve the coefficients table.

fit_summary$coefficients
#>                                  Estimate Std. Error       df     t value
#> (Intercept)                   30.96769899  0.8293349 187.9132 37.34040185
#> RACEBlack or African American  1.50464863  0.6206596 169.9454  2.42427360
#> RACEWhite                      5.61309565  0.6630909 158.8700  8.46504747
#> ARMCDTRT                       3.77555734  1.0762774 146.2690  3.50797778
#> AVISITVIS2                     4.82858803  0.8017144 143.6593  6.02282805
#> AVISITVIS3                    10.33317002  0.8224414 155.6572 12.56401918
#> AVISITVIS4                    15.05255715  1.3128602 138.3916 11.46546844
#> ARMCDTRT:AVISITVIS2           -0.01737409  1.1291645 138.3926 -0.01538668
#> ARMCDTRT:AVISITVIS3           -0.66753189  1.1865359 158.2106 -0.56258887
#> ARMCDTRT:AVISITVIS4            0.63094392  1.8507884 129.6377  0.34090549
#>                                   Pr(>|t|)
#> (Intercept)                   7.122405e-89
#> RACEBlack or African American 1.638725e-02
#> RACEWhite                     1.605553e-14
#> ARMCDTRT                      6.001485e-04
#> AVISITVIS2                    1.366921e-08
#> AVISITVIS3                    1.927523e-25
#> AVISITVIS4                    8.242710e-22
#> ARMCDTRT:AVISITVIS2           9.877459e-01
#> ARMCDTRT:AVISITVIS3           5.745112e-01
#> ARMCDTRT:AVISITVIS4           7.337266e-01

Other model parameters and metadata available in the summary object is as follows:

str(fit_summary)
#> List of 14
#>  $ cov_type        : chr "us"
#>  $ reml            : logi TRUE
#>  $ n_groups        : int 1
#>  $ n_theta         : int 10
#>  $ n_subjects      : int 197
#>  $ n_timepoints    : int 4
#>  $ n_obs           : int 537
#>  $ beta_vcov       : num [1:10, 1:10] 0.688 -0.207 -0.163 -0.569 -0.422 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:10] "(Intercept)" "RACEBlack or African American" "RACEWhite" "ARMCDTRT" ...
#>   .. ..$ : chr [1:10] "(Intercept)" "RACEBlack or African American" "RACEWhite" "ARMCDTRT" ...
#>  $ varcor          : num [1:4, 1:4] 40.73 14.27 5.14 13.53 14.27 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:4] "VIS1" "VIS2" "VIS3" "VIS4"
#>   .. ..$ : chr [1:4] "VIS1" "VIS2" "VIS3" "VIS4"
#>  $ method          : chr "Satterthwaite"
#>  $ coefficients    : num [1:10, 1:5] 30.97 1.5 5.61 3.78 4.83 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:10] "(Intercept)" "RACEBlack or African American" "RACEWhite" "ARMCDTRT" ...
#>   .. ..$ : chr [1:5] "Estimate" "Std. Error" "df" "t value" ...
#>  $ n_singular_coefs: int 0
#>  $ aic_list        :List of 4
#>   ..$ AIC     : num 3407
#>   ..$ BIC     : num 3440
#>   ..$ logLik  : num -1694
#>   ..$ deviance: num 3387
#>  $ call            : language fit_mmrm(formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),      data = "fev_data", weights = weights| __truncated__
#>  - attr(*, "class")= chr "summary.mmrm"

Other components

Specific model quantities not supported by methods can be retrieved with the component() function. The default will output all supported components.

For example, a user may want information about convergence:

component(fit, name = c("convergence", "evaluations", "conv_message"))
#> $convergence
#> [1] 0
#> 
#> $evaluations
#> function gradient 
#>       17       17 
#> 
#> $conv_message
#> [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

or the original low-level call:

component(fit, name = "call")
#> fit_mmrm(formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | 
#>     USUBJID), data = "fev_data", weights = weights, reml = reml, 
#>     control = control)

the user could also ask for all provided components by not specifying the name argument.

component(fit)

Lower level functions

Low-level mmrm

The lower level function which is called by mmrm() is fit_mmrm(). This function is exported and can be used directly. It is similar to mmrm() but lacks some post-processing and support for Satterthwaite and Kenward-Roger calculations. It may be useful for other packages that want to fit the model without the adjustment calculations.

fit_mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data,
  weights = rep(1, nrow(fev_data)),
  reml = TRUE,
  control = mmrm_control()
)
#> mmrm fit
#> 
#> Formula:     FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID)
#> Data:        fev_data (used 537 observations from 197 subjects with maximum 4 
#> timepoints)
#> Covariance:  unstructured (10 variance parameters)
#> Method:      REML
#> Deviance:    3387.373
#> 
#> Coefficients: 
#>                   (Intercept) RACEBlack or African American 
#>                   30.96769899                    1.50464863 
#>                     RACEWhite                      ARMCDTRT 
#>                    5.61309565                    3.77555734 
#>                    AVISITVIS2                    AVISITVIS3 
#>                    4.82858803                   10.33317002 
#>                    AVISITVIS4           ARMCDTRT:AVISITVIS2 
#>                   15.05255715                   -0.01737409 
#>           ARMCDTRT:AVISITVIS3           ARMCDTRT:AVISITVIS4 
#>                   -0.66753189                    0.63094392 
#> 
#> Model Inference Optimization:
#> Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmch

Hypothesis testing

This package supports estimation of one- and multi-dimensional contrasts (t-test and F-test calculation) with the df_1d() and df_md() functions. Both functions utilize the chosen adjustment method from the initial mmrm call for the calculation of degrees of freedom and (for Kenward-Roger methods) the variance estimates for the test-statistics.

One-dimensional contrasts

Compute the test of a one-dimensional (vector) contrast for a mmrm object with Satterthwaite degrees of freedom.

fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)

contrast <- numeric(length(component(fit, "beta_est")))
contrast[3] <- 1

df_1d(fit, contrast)
#> $est
#> [1] 5.643565
#> 
#> $se
#> [1] 0.6656093
#> 
#> $df
#> [1] 157.1382
#> 
#> $t_stat
#> [1] 8.478795
#> 
#> $p_val
#> [1] 1.564869e-14

This works similarly when choosing a Kenward-Roger adjustment:

fit_kr <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data,
  method = "Kenward-Roger"
)

df_1d(fit_kr, contrast)
#> $est
#> [1] 5.643565
#> 
#> $se
#> [1] 0.6740941
#> 
#> $df
#> [1] 157.1382
#> 
#> $t_stat
#> [1] 8.372073
#> 
#> $p_val
#> [1] 2.931654e-14

We see that because this is a one-dimensional contrast, the degrees of freedoms are identical for Satterthwaite and Kenward-Roger. However, the standard errors are different and therefore the p-values are different.

Multi-dimensional contrasts

Compute the test of a multi-dimensional (matrix) contrast for the above defined mmrm object with Satterthwaite degrees of freedom:

contrast <- matrix(data = 0, nrow = 2, ncol = length(component(fit, "beta_est")))
contrast[1, 2] <- contrast[2, 3] <- 1

df_md(fit, contrast)
#> $num_df
#> [1] 2
#> 
#> $denom_df
#> [1] 165.5553
#> 
#> $f_stat
#> [1] 36.91143
#> 
#> $p_val
#> [1] 5.544575e-14

And for the Kenward-Roger adjustment:

df_md(fit_kr, contrast)
#> $num_df
#> [1] 2
#> 
#> $denom_df
#> [1] 165.5728
#> 
#> $f_stat
#> [1] 35.99422
#> 
#> $p_val
#> [1] 1.04762e-13

We see that for the multi-dimensional contrast we get slightly different denominator degrees of freedom for the two adjustment methods.

Support for emmeans

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.

fit <- mmrm(
  formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)

if (require(emmeans)) {
  emmeans(fit, ~ ARMCD | AVISIT)
}
#> Loading required package: emmeans
#> AVISIT = VIS1:
#>  ARMCD emmean    SE  df lower.CL upper.CL
#>  PBO     33.3 0.757 149     31.8     34.8
#>  TRT     37.1 0.764 144     35.6     38.6
#> 
#> AVISIT = VIS2:
#>  ARMCD emmean    SE  df lower.CL upper.CL
#>  PBO     38.2 0.608 150     37.0     39.4
#>  TRT     41.9 0.598 146     40.7     43.1
#> 
#> AVISIT = VIS3:
#>  ARMCD emmean    SE  df lower.CL upper.CL
#>  PBO     43.7 0.462 131     42.8     44.6
#>  TRT     46.8 0.507 130     45.8     47.8
#> 
#> AVISIT = VIS4:
#>  ARMCD emmean    SE  df lower.CL upper.CL
#>  PBO     48.4 1.189 134     46.0     50.7
#>  TRT     52.8 1.188 133     50.4     55.1
#> 
#> Results are averaged over the levels of: RACE 
#> Confidence level used: 0.95

Acknowledgments

The mmrm package is based on previous work internal in Roche, namely the tern and tern.mmrm packages which were based on lme4. The work done in the rbmi package has been important since it used glmmTMB for fitting MMRMs.

We would like to thank Ben Bolker from the glmmTMB team for multiple discussions when we tried to get the Satterthwaite degrees of freedom implemented with glmmTMB (see https://github.com/glmmTMB/glmmTMB/blob/satterthwaite_df/glmmTMB/vignettes/satterthwaite_unstructured_example2.Rmd). Also Ben helped us significantly with an example showing how to use TMB for a random effect vector (https://github.com/bbolker/tmb-case-studies/tree/master/vectorMixed).