In longitudinal clinical trials the number and the collection times of the repeated measurements are oftentimes the same for all participants. Consequently, there is usually enough replications at each measurement occasion and within each treatment group, in the case of moderately sized or large randomized trials, to estimate the true underlying covariance matrix. Mixed-model repeated measures (MMRM) – a likelihood-based regression model for modeling the mean with time included categorically – has been generally recommended for the primary analysis of continuous outcomes with longitudinal trial designs with the use of unstructured (UN) modeling of within-patient covariances [1]. With the focus on modeling the mean response with time included categorically, the corresponding generalized estimating equations (GEE) method can provide an alternative (i.e., marginal models with an identity link). When the longitudinal measurements are multivariate Gaussian, the solution to the estimating equations is the generalized least squares estimator of the regression coefficients from MMRM [2]. Like MMRM, GEE can handle missing data easily, albeit provided that the longitudinal outcomes satisfy the missing completely at random (MCAR) assumption [3].

A compelling argument can sometimes be made to support the MCAR assumption in well-conducted trials in rare diseases. Patient centeredness is central to trial design. Involving patients, patient advocates, and family members in trial design can increase retention and ensure adequate assessment of trial outcomes [4, 5]. Compared with non-rare disease trials, trials in rare diseases tend to enroll fewer patients, address a high unmet medical need, and use patient-centered or surrogate endpoints to assess efficacy [6, 7]. Such design characteristics can mitigate dropouts and incomplete data. The use of external controls, shorter placebo-phase duration, and designs that assure participants that they will receive the experimental treatment favor patient participation and retention and further reduce the risk of dropout [8, 9]. However, if MCAR is not clinically defensible and dropout is assumed to be missing at random (MAR), the estimating equations need to be weighted by the propensity to dropout [10]. Whereas it is generally appreciated that MMRM provides a valid framework under the MAR assumption, it is not well recognized that this assumption does not hold if the covariance pattern has not been correctly specified [11]; and, therefore, may lead to a false sense of security with its use.

Little priori information about the covariance pattern is available in rare diseases. Furthermore, there are no reliable ways to identify the best covariance model in small samples [12]. If the data allow it, an UN covariance pattern should be considered; it minimizes the risk of a misspecified covariance structure. If convergence problems in numerical optimization are encountered, strategies to achieve convergence should be considered before adopting a parsimonious covariance pattern [1, 13, 14]. In some controlled or single-arm trials in rare diseases, however, there simply are not enough participants relative to the number of time points to estimate the covariance terms with an UN covariance pattern. In either case, and because we regard the correlation structure as a nuisance, the compound symmetry (CS) covariance pattern could pragmatically be considered in trials in rare diseases, which recognizes the dependence at the cost of only one extra covariance parameter.

Modeling time categorically imposes no structure on the mean response over time. It is unlikely, therefore, that the model for the mean response is misspecified; MMRM and GEE estimates of the regression coefficients will be valid. However, because the CS covariance pattern will likely be misspecified, the model-based standard errors (SEs) for MMRM and GEE estimates of the regression coefficients will be invalid [11, 15]. This is true in small and large samples [16]. With both methods, valid SEs of the regression coefficients can be based on the sandwich variance estimator, which are robust to any misspecification of the covariance [15]. Of course, the sandwich variance estimator is based on large-sample properties and its performance deteriorates for small samples [17–19]. Small-sample corrections for these SEs are needed to correct for bias.

Researchers have developed and compared a variety of bias-corrected SEs across a range of simulated scenarios, albeit separately for MMRM and GEE [18–27]. Furthermore, as opposed to relying on large-sampling properties of the regression coefficients when building confidence intervals (CIs) and testing hypotheses for the regression coefficients, t − and F-distributions using degrees of freedom (DOF) approximations to account for uncertainty in estimation of the regression coefficients can be used. Because it is not always clear what defines a small sample, a suitable bias-corrected strategy should be considered for trials conducted in rare diseases along with the choice to use either MMRM or GEE.

While there are numerous statistical publications on using these methods spanning decades, there is a dearth of literature in the rare disease community that uses real data to collectively report on both methods applied to continuous outcomes collected longitudinally from trials in rare diseases that is succinctly presented in an accessible manner. Analysts in rare disease research may be confused on which method to apply in practice and how to use the chosen method correctly. Coupled with the need to specify the primary analysis method in a protocol before the trial begins [28], a case studies approach to prespecifying a modeling strategy in rare diseases was taken. Herein, we considered three unique, albeit representative longitudinal trial designs with continuous efficacy outcomes conducted in rare diseases. We report the degree of missing data, evaluate the missing data pattern, and assess common characteristics of the empirical covariance pattern of the repeated measurements. We then simulated data directly from these trials to share with other analysts in the rare disease community a concise review of the operating characteristics of these methods based on real data. We provide a succinct summary and practical advice to complement the statistical literature on when and how to apply these methods in practice.

Within each trial, Table 4 shows the proportion of missing data, results from applying Little’s MCAR Test, and the empirical covariance pattern of the repeated measurements. Two participants in UX001-CL301 and one participant in trial UX023-CL304 discontinued the trial prematurely, while none of the participants in UX023-CL301 discontinued early. The reasons for study discontinuation were subject non-compliance (both participants in UX001-CL301) and withdrawal of consent (single participant in UX023-CL304). In total, 1.4%, 0.5%, and 1.2% of the measurement occasions in UX001-CL301, UX023-CL301, and UX023-CL304, respectively, were missing outcome data. In each trial, no missing outcome data pattern occurred in >1 participant such that the missing outcome data resembled a random sample of all the outcome data. Coupled with the clinical plausibility of MCAR and the results from applying Little’s MCAR test (all P>0.19), the longitudinal outcome data within each trial reasonably satisfied the MCAR assumption. Based on the empirical covariances for each outcome, the variances and the pairwise covariances were not constant over time.

Appendix Table 2A shows the operating characteristics associated with each method, with and without the sandwich variance estimator and small-sample correction, based on 100,000 simulated datasets each from UX001-CL301 and UX023-CL301 and 500,000 from UX023-CL304. The prespecified GEE method assumed a CS covariance pattern with the use of the sandwich variance estimator without a small-sample correction and inference based on asymptotic normality, while the data were generated from MMRM with UN covariance pattern. On average, in any given dataset simulated from UX001-CL301, UX023-CL301, and UX023-CL304, 1.4%, 0.5%, and 1.2% of the longitudinal measurements were MCAR.

MMRM, a response profile analysis, arguably dates back to Greenhouse and Geisser (1959) [38], and has been widely adopted as the primary analysis for analyzing repeatedly measured continuous outcomes in randomized clinical trials (RCTs). For an example of a recently completed phase 3 trial (N = 65; NCT03399786) in a rare disease that applied MMRM as the primary analysis to assess a treatment effect at a single time point, albeit without the sandwich variance estimator and small-sample correction, see Raal et al (2020) [39]. Since the two initial companion papers on GEE for longitudinal data [3, 40] coupled with the connections drawn between GEE and likelihood-based methods [41] – and the widely available software notwithstanding – GEE has infrequently been used as the primary analysis for continuous outcomes in RCTs. Based on a recent comprehensive evaluation of statistical methods applied as the primary analysis to repeatedly measured continuous outcomes from RCTs published in 2019, Ren et al (2022) [42] estimated 4% (4/96) adopted GEE to assess a treatment effect at a single time point compared with 36% (35/96) that used MMRM; the remaining 59% (57/96) oversimplified the longitudinal design and applied conventional methods such as t-tests and analysis of covariance to assess a treatment effect at a single time point. One rare disease trial (N = 53; NCT02208687) that adopted in primary analysis GEE with the sandwich variance estimator was designed to investigate the effectiveness at a single time point of a self-management group program to improve social participation and endurance in patients with neuromuscular disease with chronic fatigue [43]. Either method can be considered for modeling longitudinal changes in the mean response in trials conducted in rare diseases with a parsimonious CS covariance structure; however, care needs to be exercised with both methods to ensure its correct use.

Inferences about longitudinal changes in mean response and its relationship to treatment arm are sensitive to the chosen covariance model. This is true for both MMRM and GEE when modeling the mean with time included categorically, and despite the covariances being treated as nuisance parameters with GEE. In rare disease trials, it is common, however, to pragmatically prespecify a parsimonious CS covariance pattern due to the smaller sample sizes or to prespecify a decision to adopt a CS covariance pattern if an UN covariance pattern was inestimable. Based on our case studies, the variances and covariances of the longitudinal outcomes were not constant over time. When the CS covariance pattern does not closely approximate the true but unknown covariance pattern, incorrect SEs of the regression coefficients will result. As seen with our case studies, to obtain valid SEs, the sandwich variance estimator coupled with a small-sample correction to correct for bias, must be considered with both methods.

With GEE, the sandwich variance estimator – as opposed to model-based SEs – is the default setting in most statistical packages, however, the opposite is true when fitting MMRM. Consider, for example, SAS and R. In SAS, the user needs to specify the EMPIRICAL option in the PROC GLIMMIX (or MIXED) statement to apply the sandwich variance estimator when fitting MMRM, while the sandwich variance estimator is the default when using GENMOD to fit the corresponding GEE. When fitting MMRM in R using the gls function in the nlme package [44], the analyst would inconveniently appeal to the clubSandwich package [45] and the function vcovCR, which returns a sandwich estimate of the variance-covariance matrix of the regression coefficients. The recently developed mmrm package [13] currently does not have an option to apply sandwich variance estimators. To fit the corresponding GEE in R, the user can appeal to the gee [46] or geepack [47] packages, and in both packages the sandwich variance estimator is the default. Defaulting to the sandwich variance estimator with GEE is not surprising. The properties of the sandwich variance estimator has been well studied with GEE [19, 48–51]. If there is misspecification of the covariance, however, not using the sandwich variance estimator can also result in misleading inferences for the regression coefficients in MMRM as seen in our case studies.

Some small-sample corrections to the sandwich variance estimators are available in SAS and R to correct for bias along with some DOF approximations to the t distribution. Pustejovsky and Tipton (2017) [24] and Wang and coauthors (2016) [25] provide a comprehensive review of small-sample corrections available in R and guidance on which approach to use for the methods considered in this article. Gosho et al (2021) [22] provide a practical review of the small-sample corrections for MMRM that can be obtained using the GLIMMIX procedure directly or programmatically using IML.

Longitudinal responses will invariably be incomplete. Unlike MCAR, when the missing response data mechanism is MAR, subjects with missing response data will no longer be a random subset of the sample. It is well known that MMRM provides a valid framework under MAR, a more defensible assumption than MCAR in any given trial, albeit provided that the covariance pattern has been correctly specified [1]. If estimating an UN covariance pattern is prohibitive due to the small sample, resulting in misspecification, inferences about the mean response over time will be invalid when dropout is MAR. This point further underscores the importance of considering the sandwich variance estimator with MMRM (± a small-sample correction). When dropout is MAR as opposed to MCAR, the validity of the analysis with the standard GEE method based on available observations will be compromised [14]. Inverse probability weighted (IPW) [10] GEE should be adopted to adjust the analysis for the propensity to dropout; otherwise, GEE will yield biased estimates of the coefficients in the model part of the mean response. The sandwich variance estimator when using IPW methods can be constructed to adjust for estimation of the weights. Because most statistical procedures allow for the inclusion of sampling weights, it is straightforward to apply IPW-GEE. While unlikely in a well-controlled longitudinal clinical trial [1], if the longitudinal responses are missing not at random (MNAR), none of the methods considered in this article will be valid [11].

In our case studies, the MCAR assumption was clinically defensible. The number of measurement occasions including baseline was ≤7, few patients discontinued prematurely, no missing outcome data pattern occurred in >1 participant, and the degree of missing outcome data was <1.5%. While the statistical test of the MCAR assumption likely had low power [52], its use combined with the empirical evidence provided justification to use GEE. Due to the hierarchy of dropout mechanisms [53], MMRM was also justified.

There were research limitations. First, we only considered GLIMMIX to estimate MMRM and the SAS/IML GEE macro, which reproduces the method-of-moments approach to estimation with GENMOD. SAS is widely used software and the primary analysis in each trial used GENMOD. Second, there are numerous bias-corrections to sandwich variance estimators in the literature but only a few are currently available with GLIMMIX and none with GENMOD. The choice of which bias-correction to use would depend on several factors [19, 24, 25] and should be prespecified. In this article, we only considered the bias-corrected SEs recommended by Mancl and DeRouen (2001) [19]. Third, for t-based inference, we only considered the between-within method DOF approximation for small-sample inference recommended by Schluchter and Elashoff (1990) [36]. When empirical or bias-corrected SEs are used with GLIMMIX, DOF approximation methods are limited to between-within, containment, and residual, where the latter two are not small-sample approximations; furthermore, GENMOD only performs inference assuming asymptotic normality. The determination of which DOF approximation to use is an active area of research, with many methods not accessible in standard SAS and R routines. Lastly, the research was limited to three case studies and, in each, <1.5% of the measurement occasions were missing outcome data. Although these case studies represented typical longitudinal designs in rare diseases, it was not an exhaustive representation and the conclusions may not necessarily extend to other design variations. In particular, dropout other than MCAR was not considered in the simulation study. However, to understand whether the findings held in the presence of moderate missingness defined as 10% of the measurement occasions being MCAR, the simulation study was repeated. The simulation-based conclusions when using a CS covariance pattern remained; however, for the smallest trial (UX023-CL304; N = 14), estimation was unstable or not always possible when assuming an unstructured covariance pattern. In any case, the intent of the current research was not to pressure test these methods; rather to elucidate the issues and the steps to address the issues using real datasets in rare diseases to aid an analyst working in a rare disease in selecting an appropriate modeling strategy.

Based on a review of three case studies recently conducted in rare diseases, the MCAR assumption was plausible and the missingness low. When modeling the mean response with time as a categorical variable, MMRM with UN covariance coupled with the Kenward-Roger standard error adjustment should be used if the data allow it. If not, then a parsimonious CS covariance structure can be considered. In the two case studies that exhibited nonconstant variance/covariances over time, the sandwich variance estimator and small-sample correction with t-based inference using the between-within DOF should be considered with both MMRM and GEE. If the CS pattern was a good approximation, as seen in the other case study, then model-based standard errors with t-based inference using the between-within DOF performed well with both methods.

At the time the research was conducted, David Zahrieh, Yi Wang, and Tony Koutsoukos were employees at Ultragenyx Pharmaceutical Inc. This research reflects the views of the authors and should not be construed to represent Ultragenyx’s views or policies.

Data sharing is not applicable to this article as no new data were created or analyzed in this research.

This work was presented as a Poster at the 2023 Joint Statistical Meetings on August 7, 2023.