From the high-level mmrm()
interface, common changes to the default function call can be
specified.
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.
The starting values will affect the optimization result. A better
starting value usually can make the optimization more efficient. In
mmrm we provide two starting value functions, one is
std_start and the other is emp_start.
std_start will try to use the identity matrix as the
covariance, however there are convergence problems for ar1
and ar1h if the identity matrix is provided, thus for these
two covariance structures we use \(\rho=0.5\) instead. emp_start
will try to use the empirical covariance matrix of the residuals of the
ordinary least squares model as the starting value for unstructured
covariance structure. If some timepoints are missing from data, identity
matrix will be used for that submatrix. The correlation between existing
and non-existing timepoints are set to 0.
As the starting values will affect the result, please be cautious on choosing the starting values.
Here we provide an example where the std_start fails. In
the following code chunk, we will create a dummy dataset for mmrm
analysis.
gen_data <- function(
n = 100,
mu = -100 / 52,
delta = 50 / 52,
mua = 2000,
sigmaa = 300,
sigmab = 60,
corab = 0.2,
sigma = 10,
times = c(0, 2, 6, 12, 24, 36, 52, 70, 88, 104)) {
nt <- length(times)
out <- data.frame(
pts = rep(seq_len(n * 2), each = nt),
trt = rep(c("Treatment", "Placebo"), rep(n * nt, 2)),
time = rep(times, n * 2)
)
covab <- corab * sigmaa * sigmab # cov between a and b
cov <- matrix(c(sigmaa^2, covab, covab, sigmab^2), ncol = 2) # Cov matrix for the slope and intercept
si <- rbind(
MASS::mvrnorm(n, mu = c(mua, mu + delta), Sigma = cov),
MASS::mvrnorm(n, mu = c(mua, mu + delta), Sigma = cov)
)
idx <- rep(seq_len(n * 2), each = nt)
out$fev1 <- si[idx, 1] + si[idx, 2] * times + rnorm(n * nt * 2, sd = sigma)
out$trt <- factor(out$trt)
out$time <- factor(out$time)
out$pts <- factor(out$pts)
return(out)
}
set.seed(123)
out <- gen_data()In the generated data, the variance is not in the same scale across visits.
vapply(split(out$fev1, out$time), sd, FUN.VALUE = 1)## 0 2 6 12 24 36 52 70
## 278.6079 319.0589 482.4172 799.9107 1491.1440 2194.5776 3140.0768 4204.9355
## 88 104
## 5272.6041 6221.2195Using emp_start as the starting value, mmrm
will converge fast.
mmrm(fev1 ~ trt * time + us(time | pts), data = out, start = emp_start)## mmrm fit
##
## Formula: fev1 ~ trt * time + us(time | pts)
## Data: out (used 2000 observations from 200 subjects with maximum 10
## timepoints)
## Covariance: unstructured (55 variance parameters)
## Inference: REML
## Deviance: 19154.63
##
## Coefficients:
## (Intercept) trtTreatment time2
## 1962.6980059 11.3831958 0.2130684
## time6 time12 time24
## -1.6901336 -1.1984278 -7.3954489
## time36 time52 time70
## -11.4078896 -15.8040921 -22.5556525
## time88 time104 trtTreatment:time2
## -28.2068896 -33.6608068 6.0436315
## trtTreatment:time6 trtTreatment:time12 trtTreatment:time24
## 27.3686373 49.4246567 107.3638488
## trtTreatment:time36 trtTreatment:time52 trtTreatment:time70
## 161.4310444 233.6438342 316.7387101
## trtTreatment:time88 trtTreatment:time104
## 397.9895966 471.8871913
##
## Model Inference Optimization:
## Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmchHowever, if we use std_start, there will be convergence
problems. We can also force a specific optimization algorithm and add
control details, here e.g. choosing nlminb with increased
maximum number of function evaluations and iterations.
mmrm(
fev1 ~ trt * time + us(time | pts),
data = out,
start = std_start,
optimizer = "nlminb",
optimizer_control = list(eval.max = 1000, iter.max = 1000)
)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)
## Inference: 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*epsmchUsers 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,
nlminb and other user-defined custom optimizers.
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)
## Inference: REML
## Deviance: 3387.373
##
## Coefficients:
## (Intercept) RACEBlack or African American
## 30.9676894 1.5046748
## RACEWhite ARMCDTRT
## 5.6131058 3.7755440
## AVISITVIS2 AVISITVIS3
## 4.8285859 10.3331765
## AVISITVIS4 ARMCDTRT:AVISITVIS2
## 15.0525711 -0.0173539
## ARMCDTRT:AVISITVIS3 ARMCDTRT:AVISITVIS4
## -0.6675207 0.6309568
##
## Model Inference Optimization:
## Converged with code 0 and message: No message provided.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)
## Inference: ML
## Deviance: 3536.989
##
## Coefficients:
## (Intercept) RACEBlack or African American
## 31.4207078 0.5357236
## RACEWhite ARMCDTRT
## 5.4546329 3.4305211
## 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*epsmchThe 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 coded as a factor.
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)
## Inference: 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*epsmchGrouped 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.822313 3.600996 3.600996 3.600996
## VIS2 3.600996 37.822313 3.600996 3.600996
## VIS3 3.600996 3.600996 37.822313 3.600996
## VIS4 3.600996 3.600996 3.600996 37.822313
##
## $TRT
## VIS1 VIS2 VIS3 VIS4
## VIS1 49.58031 10.98097 10.98097 10.98097
## VIS2 10.98097 49.58031 10.98097 10.98097
## VIS3 10.98097 10.98097 49.58031 10.98097
## VIS4 10.98097 10.98097 10.98097 49.58031We can see that the estimated covariance matrices are different in
different ARMCD groups.
In additional to the residual and Between-Within degrees of freedom,
both Satterthwaite and Kenward-Roger adjustment methods are available.
The default is Satterthwaite adjustment of the degrees of freedom. To
use e.g. the Kenward-Roger adjustment of the degrees of freedom as well
as the coefficients covariance matrix, use the method
argument:
A list of all allowed method is
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
## Vcov 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.6006For 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",
vcov = "Kenward-Roger-Linear"
)This is due to the different parametrization of the unstructured covariance matrix, see the Kenward-Roger vignette for details.
There are multiple variance-covariance estimator available for the coefficients, including:
Please note that, not all combinations of variance-covariance for coefficients and method of degrees of freedom are possible, e.g. “Kenward-Roger” and “Kenward-Roger-Linear” are available only when the degrees of freedom method is “Kenward-Roger”.
Details can be found in Coefficients Covariance Matrix Adjustment vignette and Weighted Least Square Empirical Covariance.
An example of using other variance-covariance is:
fit_emp <- mmrm(
formula = FEV1 ~ RACE + ARMCD * AVISIT + us(AVISIT | USUBJID),
data = fev_data,
method = "Satterthwaite",
vcov = "Empirical"
)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.
## Additional attributes including contrasts are lost.
## To avoid this behavior, make sure use `drop_visit_levels = FALSE`.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.00000000cov2cor(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.00000000We 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.
In some scenarios, we may want to control how the design matrix is
constructed in finer detail. This can be especially of interest if there
are factor levels not present in the current dataset, which nonetheless
are valid levels that should not result in later errors (for example,
when using predict functions). In the example below,
mmrm might not run if only one factor level is
available.
ex_data <- fev_data[!(fev_data$RACE == "White" & fev_data$SEX == "Male"),]
fit_subgroup1 <- mmrm(
formula = FEV1 ~ SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
data = ex_data[ex_data$RACE != "White",],
vcov = "Asymptotic"
)
fit_subgroup2 <- mmrm(
formula = FEV1 ~ SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
data = ex_data[ex_data$RACE == "White",],
vcov = "Asymptotic"
)## Some factor levels are dropped due to singular design matrix: SEX## Error in `contrasts<-`:
## ! contrasts can be applied only to factors with 2 or more levelsemmeans::emmeans(fit_subgroup1, ~ARMCD|AVISIT, data=fev_data, weights="proportional")## mmrm() registered as emmeans extension## AVISIT = VIS1:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 31.8 0.812 106.9 30.2 33.4
## TRT 35.2 0.916 106.1 33.4 37.0
##
## AVISIT = VIS2:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 36.9 0.671 106.8 35.5 38.2
## TRT 40.4 0.720 104.2 38.9 41.8
##
## AVISIT = VIS3:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 41.9 0.515 94.6 40.9 43.0
## TRT 44.9 0.631 95.2 43.6 46.1
##
## AVISIT = VIS4:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 46.3 1.340 96.0 43.6 48.9
## TRT 50.9 1.470 95.4 47.9 53.8
##
## Results are averaged over the levels of: SEX
## Confidence level used: 0.95emmeans::emmeans(fit_subgroup2, ~ARMCD|AVISIT, data=fev_data, weights="proportional")## Error:
## ! object 'fit_subgroup2' not foundBy setting the contrasts argument specifically, just like in
lm(), we can avoid this error:
contr_mat <- contr.sum(nlevels(ex_data$SEX))
rownames(contr_mat) <- levels(fev_data$SEX)
contr_mat## [,1]
## Male 1
## Female -1Now the subgroup analysis code will run correctly, even when
emmeans averages over the baseline values that aren’t
subgroup-specific:
fit_subgroup1 <- mmrm(
formula = FEV1 ~ SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
data = ex_data[ex_data$RACE == "White",],
vcov = "Asymptotic",
contrasts = list(SEX = contr_mat)
)
fit_subgroup2 <- mmrm(
formula = FEV1 ~ SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
data = ex_data[ex_data$RACE != "White",],
vcov = "Asymptotic",
contrasts = list(SEX = contr_mat)
)
emmeans::emmeans(fit_subgroup1, ~ARMCD|AVISIT, data=fev_data, weights="proportional")## AVISIT = VIS1:
## ARMCD emmean SE df asymp.LCL asymp.UCL
## PBO nonEst NA NA NA NA
## TRT nonEst NA NA NA NA
##
## AVISIT = VIS2:
## ARMCD emmean SE df asymp.LCL asymp.UCL
## PBO nonEst NA NA NA NA
## TRT nonEst NA NA NA NA
##
## AVISIT = VIS3:
## ARMCD emmean SE df asymp.LCL asymp.UCL
## PBO nonEst NA NA NA NA
## TRT nonEst NA NA NA NA
##
## AVISIT = VIS4:
## ARMCD emmean SE df asymp.LCL asymp.UCL
## PBO nonEst NA NA NA NA
## TRT nonEst NA NA NA NA
##
## Results are averaged over the levels of: SEX
## Confidence level used: 0.95emmeans::emmeans(fit_subgroup2, ~ARMCD|AVISIT, data=fev_data, weights="proportional")## AVISIT = VIS1:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 31.8 0.812 106.9 30.2 33.4
## TRT 35.2 0.916 106.1 33.4 37.0
##
## AVISIT = VIS2:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 36.9 0.671 106.8 35.5 38.2
## TRT 40.4 0.720 104.2 38.9 41.8
##
## AVISIT = VIS3:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 41.9 0.515 94.6 40.9 43.0
## TRT 44.9 0.631 95.2 43.6 46.1
##
## AVISIT = VIS4:
## ARMCD emmean SE df lower.CL upper.CL
## PBO 46.3 1.340 96.0 43.6 48.9
## TRT 50.9 1.470 95.4 47.9 53.8
##
## Results are averaged over the levels of: SEX
## Confidence level used: 0.95