Benchmarking
Next, the MMRM fitting procedures are compared using the FEV and BCVAdatasets. FEV1 measurements are modeled as a function of race, treatmentarm, visit number, and the interaction between the treatment arm and thevisit number. Change in BCVA is assumed to be a function of race,baseline BCVA, treatment arm, visit number, and the treatment–visitinteraction. In both datasets, repeated measures are modeled using anunstructured covariance matrix. The implementations’ convergence timesare evaluated first, followed by a comparison of their estimates.Finally, we fit these procedures on simulated BCVA-like data to assessthe impact of missingness on convergence rates.
Convergence Times
FEV Data
Themmrm,PROC GLIMMIX,gls,lmer, andglmmTMB functions are applied to theFEV dataset 10 times. The convergence times are recorded for eachreplicate and are reported in the table below.
| Implementation | Median | First Quartile | Third Quartile |
|---|---|---|---|
| mmrm | 56.15 | 55.76 | 56.30 |
| PROC GLIMMIX | 100.00 | 100.00 | 100.00 |
| lmer | 247.02 | 245.25 | 257.46 |
| gls | 687.63 | 683.50 | 692.45 |
| glmmTMB | 715.90 | 708.70 | 721.57 |
It is clear from these results thatmmrm convergessignificantly faster than other R functions. Though not demonstratedhere, this is generally true regardless of the sample size andcovariance structure used.mmrm is faster thanPROC GLIMMIX.
BCVA Data
The MMRM implementations are now applied to the BCVA dataset 10times. The convergence times are presented below.
| Implementation | Median | First Quartile | Third Quartile |
|---|---|---|---|
| mmrm | 3.36 | 3.32 | 3.46 |
| glmmTMB | 18.65 | 18.14 | 18.87 |
| PROC GLIMMIX | 36.25 | 36.17 | 36.29 |
| gls | 164.36 | 158.61 | 165.93 |
| lmer | 165.26 | 157.46 | 166.42 |
We again find thatmmrm produces the fastest convergencetimes on average.
Marginal Treatment Effect Estimates Comparison
We next estimate the marginal mean treatment effects for each visitin the FEV and BCVA datasets using the MMRM fitting procedures. All Rimplementations’ estimates are reported relative toPROC GLIMMIX’s estimates. Convergence status is alsoreported.
FEV Data
The R procedures’ estimates are very similar to those output byPROC GLIMMIX, thoughmmrm andglsgenerate the estimates that are closest to those produced when usingSAS. All methods converge using their default optimizationarguments.
BCVA Data
mmrm,gls andlmer produceestimates that are virtually identical toPROC GLIMMIX’s,whileglmmTMB does not. This is likely explained byglmmTMB’s failure to converge. Note too thatlmer fails to converge.
Impact of Missing Data on Convergence Rates
The results of the previous benchmark suggest that the amount ofpatients missing from later time points affect certain implementations’capacity to converge. We investigate this further by simulating datausing a data-generating process similar to that of the BCVA datasets,though with various rates of patient dropout.
Ten datasets of 200 patients are generated each of the followinglevels of missingness: none, mild, moderate, and high. In all scenarios,observations are missing at random. The number patients observed at eachvisit is obtained for one replicated dataset at each level ofmissingness is presented in the table below.
| none | mild | moderate | high | |
|---|---|---|---|---|
| VIS01 | 200 | 196.7 | 197.6 | 188.1 |
| VIS02 | 200 | 195.4 | 194.4 | 182.4 |
| VIS03 | 200 | 195.1 | 190.7 | 175.2 |
| VIS04 | 200 | 194.1 | 188.4 | 162.8 |
| VIS05 | 200 | 191.6 | 182.5 | 142.7 |
| VIS06 | 200 | 188.2 | 177.3 | 125.4 |
| VIS07 | 200 | 184.6 | 168.0 | 105.9 |
| VIS08 | 200 | 178.5 | 155.4 | 82.6 |
| VIS09 | 200 | 175.3 | 139.9 | 58.1 |
| VIS10 | 200 | 164.1 | 124.0 | 39.5 |
The convergence rates of all implementations for stratified bymissingness level is presented in the plot below.
mmrm,gls, andPROC GLIMMIXare resilient to missingness, only exhibiting some convergence problemsin the scenarios with the most missingness. These implementationsconverged in all the other scenarios’ replicates.glmmTMB,on the other hand, has convergence issues in the no-, mild-, andhigh-missingness datasets, with the worst convergence rate occurring inthe datasets with the most dropout. Finally,lmer isunreliable in all scenarios, suggesting that it’s convergence issuesstem from something other than the missing observations.
Note that the default optimization schemes are used for each method;these schemes can be modified to potentially improve convergencerates.
A more comprehensive simulation study using data-generating processessimilar to the one used here is outlined in thesimulations/missing-data-benchmarkssubdirectory. In addition to assessing the effect of missing data onsoftware convergence rates, we also evaluate these methods’ fit timesand empirical bias, variance, 95% coverage rates, type I error rates andtype II error rates.mmrm is found to be the most mostrobust software for fitting MMRMs in scenarios where a large proportionof patients are missing from the last time points. Additionally,mmrm has the fastest average fit times regardless of theamount of missingness. All implementations considered produce similarempirical biases, variances, 95% coverage rates, type I error rates andtype II error rates.