{'fraction_non_alphanumeric': 0.04650369528140989, 'fraction_numerical': 0.016850483229107448, 'mean_word_length': 4.40312077650817, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robustly leveraging the post-randomization information to improve precision in the analyses of randomized clinical trials

In randomized clinical trials, repeated measures of the outcome are routinely collected. The mixed model for repeated measures (MMRM) leverages the information from these repeated outcome measures, and is often used for the primary analysis to estimate the average treatment effect at the final visit. MMRM, however, can suffer from precision loss when it models the intermediate outcomes incorrectly, and hence fails to use the post-randomization information in a harmless way. In this paper, we propose a new working model, called IMMRM, that generalizes MMRM and optimizes the precision gain from covariate adjustment, stratified randomization and adjustment for intermediate outcome measures. We prove that the IMMRM estimator for the average treatment effect is consistent and asymptotically normal under arbitrary misspecification of its working model assuming missing completely at random. Under simple or stratified randomization, the IMMRM estimator is asymptotically equally or more precise than the analysis of covariance (ANCOVA) estimator and the MMRM estimator. By re-analyzing three randomized trials in the diabetes therapeutic area, we demonstrate that the IMMRM estimator has 2-24% smaller variance than ANCOVA and 5-16% smaller variance than MMRM.


Introduction
"Intermediate outcomes" in the context of clinical trials refer to the outcomes measured after treatment assignment but before the time point of interest. For example, for a trial studying the effect of a 12-month dietary plan on weight, the intermediate outcomes can be the weights at 3, 6 and 9 months after randomization. In randomized clinical trials, intermediate outcomes are routinely collected and have been used in various ways, including trial monitoring (Shih and Quan, 1999), principal stratification (Seuc et al., 2013), decision making in the interim analysis (Kunz et al., 2015), mediation analysis (Landau et al., 2018), etc. We focus on a less known purpose of intermediate outcomes, which is to robustly improve the precision of statistical inference at the time point of interest, e.g. the primary endpoint.
For the purpose of increasing precision and power in clinical trials, covariate adjustment and stratified randomization are two commonly used methods. Covariate adjustment refers to adjustment for chance imbalances in baseline variables, called covariates, among treatment groups by a regression model, and has been extensively studied and applied as a robust method to reduce variance. See Yang and Tsiatis (2001); Rubin and van der Laan (2008); Moore and van der Laan (2009a,b); Zhang (2015); Jiang et al. (2018), just to name a few, for examples. Stratified randomization (Zelen, 1974), also known as stratified permuted block randomization, is a popular randomization scheme that minimizes treatment imbalance within each prespecified randomization strata, and can also increase power (Bugni et al., 2018Wang et al., 2019b;Ye et al., 2020). According to a survey by Lin et al. (2015), stratified randomization is used by 70% of trials published in top medical journals. The recent guidance from the U.S. Food and Drug Administration (FDA, 2021) advocates using covariate adjustment and stratified randomization for improving precision. The guidance, however, pointed out that one of its limitations is to not address the use of covariate adjustment for analyzing longitudinal repeated measures data. Our goal is to adjust for intermediate outcomes on top of the covariate adjustment and stratified randomization to maximize the precision gain for the analyses of randomized clinical trials with continuous outcomes. Marschner and Becker (2001) first showed that jointly modeling short-term and long-term binary outcomes can lead to precision gain compared to modeling long-term outcomes only. Galbraith and Marschner (2003); Stallard (2010); Hampson and Jennison (2013); Zhou et al. (2018) focused on interim analyses and showed that adjustment for prognostic short-term outcomes can improve interim decisions, assuming the outcomes follow a bivariate normal distribution. All of these methods, however, rely on correctly specified parametric models.
Building on general approaches developed by (Lu and Tsiatis, 2011;van der Laan and Gruber, 2012), Qian et al. (2019) derived the non-parametric formulas for the precision gain from adjusting for baseline variables and the short-term outcomes; and Van Lancker et al. (2020) proposed a model-robust estimator that achieves such precision gain for the interim analysis.
However, their asymptotic results are limited to simple randomization, single intermediate outcome, and monotone censoring, which fail to take into account the precision gain from stratified randomization and the complexity of non-monotone censoring at multiple visits.
In this paper, we propose a new working model "Improved Mixed Model for Repeated Measures" (IMMRM), which can combine the precision gain from covariate adjustment, stratified randomization and adjustment for intermediate outcomes. The IMMRM working model is an extension of the mixed model for repeated measures (MMRM, Mallinckrod et al. 2008) to handle multiple treatment groups, treatment heterogeneity and heteroscedasticity.
To the best of our knowledge, this is the first method that fully utilizes pre-randomization variables, the randomization procedure, and post-randomization information to improve precision.
Assuming mild regularity conditions and missing completely at random, the IMMRM estimator for the average treatment effect at the time point of interest is model-robust, i.e., consistent and asymptotically normal under arbitrary misspecification of its working model. Furthermore, under simple or stratified randomization, the IMMRM estimator is asymptotically equally or more precise than the following commonly-used estimators: the analysis of covariance (ANCOVA, Yang and Tsiatis, 2001) estimator, and the MMRM estimator with or without visit-by-covariates interactions in the working model.
Our result implies that appropriately leveraging post-randomization information can lead to precision gain beyond what comes from covariate adjustment and stratified randomization.
In contrast, although MMRM involves intermediate outcomes, it can be less precise than ANCOVA, even when visit-by-covariates interactions are included in its working model. This result generalizes the finding of Schuler (2021), which showed ANCOVA outperforms MMRM under two-arm, equal, simple randomization with no missing data. In addition, we provide the necessary and sufficient condition for when intermediate outcomes provide precision gain, which is, essentially, that the intermediate outcomes provide additional information beyond covariates to the missing outcomes at the time point of interest.
In the next section, we introduce three randomized clinical trials on type 2 diabetes. In Section 3, we present our setup, notations and assumptions. In Section 4, we describe the ANCOVA estimator, MMRM estimator (with or without visit-by-covariates interactions), and our proposed IMMRM estimator. Our main result is presented in Section 5, which consists of asymptotic theory, explanations of precision gain from intermediate outcomes, and discussions on a special case of two-arm equal randomization. Simulation study and data application on three completed diabetes randomized clinical trials are given in Section 6 and 7, respectively. We provide practical recommendations and discuss future directions in Section 8.
2 Three randomized clinical trials 2.1 Trial 1: the IMAGINE-2 study The trial of "A Study in Patients With Type 2 Diabetes Mellitus (IMAGINE 2)" (NCT01435616) is a 52-week, two-arm, phase 3 randomized clinical study completed in 2014 (Davies et al., 2016). The goal of this trial was to evaluate the effect of basal insulin for treatment of type 2 diabetes in an insulin-naïve population.
Participants were randomly assigned to receive insulin peglispro (treatment, 1003 patients) or insulin glargine (control, 535 patients), with a target to achieve 2:1 randomization ratio. Randomization was stratified by baseline HbA1c (< 8.5% or ≥ 8.5%), low-density lipoprotein cholesterol (< 100mg/dL or ≥ 100mg/dL) and baseline sulphonylurea or meglitinide use (yes or not). HbA1c is a continuous measure of average blood glucose values in the prior three months. The primary outcome was the change in HbA1c at week 52 from baseline (15% missing outcomes), while intermediate outcomes were measured at week 4 (3% missing outcomes), week 12 (6% missing outcomes), week 26 (10% missing outcomes) and week 39 (13% missing outcomes). We focused on estimating the average treatment effect of the primary outcome and adjusted for intermediate outcomes, baseline HbA1c value, and randomization strata across different estimators.

Trial 2: another insulin peglispro study
The trial of "A Study of Insulin Peglispro in Participants With Type 2 Diabetes Mellitus" (NCT02106364) is a 26-week, two-arm, phase 3 randomized clinical study comparing insulin peglispro with insulin glargine for treatment of type 2 diabetes mellitus in Asian insulin-naïve population (Hirose et al., 2018). This trial studied the same product (insulin peglispro) as in Trial 1, while focused on the Asian population.
The primary outcome was the change in HbA1c at week 26 from baseline (6% missing outcomes). Intermediate outcomes at week 4 (2% missing outcomes) and week 12 (3% missing outcomes) were included to the estimation. We focused on estimating the average treatment effect of the primary outcome with adjustment for intermediate outcomes, baseline HbA1c value, and randomization strata across different estimators.
The goal of this trial was to evaluate the dose-response in efficacy and safety of four doses of tirzepatide (1mg, 5mg, 10mg and 15mg), a novel dual GIP and GLP-1 receptor agonist, in patients with type 2 diabetes compared to the placebo and dulaglutide 1.5mg (active comparator). In this study, 318 patients were equally randomized to one of the six parallel treatment groups with each group containing 51-55 patients. For the purpose of demonstration, we dropped the active comparator group from our analysis and only considered comparing four doses of tirzepatide to the placebo.
Randomization was stratified by baseline HbA1c (< 8.5% or ≥ 8.5%), metformin use (yes or no), and BMI (< 30 kg/m² or ≥ 30 kg/m²). We focused on the the body weight change from baseline, which was a secondary outcome of the study. The weight change was measured at eight post-randomization visits, which are week 1, 2, 4, 8, 12, 16, 20 and 26, with 1%, 3%, 6%, 11%, 15%, 18%, 17%, 19% missing outcomes at each visit respectively. Our data application focused on estimating the average treatment effect of the weight change at week 26 from baseline and adjusted for all intermediate outcomes, the baseline body weight, and randomization strata across different estimators.

Data generating distributions
Consider a trial where the outcome is continuous and repeatedly measured at K visits, where K is a positive integer. For each participant i = 1, . . . , n, the outcome vector at K visits is where M it = 1 if Y it is observed at visit t, and 0 otherwise. Let A i be a categorical variable taking values in {0, 1, . . . , J : J ≥ 1}, with A i = j representing that participant i is assigned to the j-th treatment group. By convention, we use A i = 0 to denote being assigned to the control group. Let X i be a vector of baseline variables with length p. Throughout, we refer to Y K as the final outcome, and (Y 1 , . . . , Y K−1 ) as the intermediate outcomes (if K > 1), for conciseness.
We use the Neyman-Rubin potential outcomes framework (Neyman et al., 1990), which is the potential outcome vector for treatment group j = 0, . . . , J. Analogously to the consistency assumption above, we assume M i = J j=0 I{A i = j}M i (j), where M i (j) is the indicator vector of whether Y i (j) would be observed at each of the K visits, if participant i was assigned to treatment group j for j = 0, . . . , J.
For each participant i, we define the complete data vector as where Y o i is the vector of observed outcomes, whose dimension may vary across participants. For example, if participant i only shows up in visits 1 and K, then . For the special case that a participant misses all post-randomization visits but still has baseline information recorded, the observed data vector is with M i being a zero vector. All estimators defined in Section 4 are functions of the observed We make the following assumptions on (W 1 , . . . , W n ): . . , n are independent, identically distributed samples from the joint distribu- 2. (Missing completely at random, MCAR) M (j) is independent of (Y (j), X) for j = 0, . . . , J, and M (0), . . . , M (J) are identically distributed.
In addition to Assumption 1, we also assume regularity conditions for the estimators defined in Section 4. As we show in the Supplementary Material, all estimators we consider (including our proposed estimators) are M-estimators (Section 5 of van der Vaart, 1998), which is defined as a zero of prespecified estimating functions. The regularity conditions are made on these estimating functions, the complete data distribution P , and the parameters involved in the estimating equations. These conditions are similar to the classical conditions given in Section 5.3 of van der Vaart (1998) for proving consistency and asymptotic normality for M-estimators under simple randomization. We provide these regularity conditions in the Supplementary Material.
The parameters of interest are the average treatment effects of the final outcome comparing each treatment group to the control group, i.e., where E is the expectation with respect to the distribution P . Our results in Section 5 also apply for estimating any linear transformation of ∆ * , e.g., the average treatment effect comparing any two treatment groups.

Simple and stratified randomization
We consider two types of treatment assignment procedures: simple randomization and stratified randomization. For j = 0, . . . , J, let π j be the target proportion of participants that are assigned to the treatment group j. We assume that J j=0 π j = 1 and π j > 0 for all treatment groups. For example, equal randomization refers to the setting where π 0 = · · · = π J .
Simple randomization allocates treatment by independent draws from a categorical distribution on A, with P (A = j) = π j for j = 0, . . . , J. Then (A 1 , . . . , A n ) are independent, identically distributed samples from this categorical distribution, and also independent of (W 1 , . . . , W n ).
Under stratified randomization, treatment assignment depends on a set of categorical baseline variables, which are called stratification variables. We use a categorical random variable S with support S = {1, . . . , R} to denote the joint levels created by all stratification variables. For example, if randomization is stratified by sex (female or male) and weight (normal, overweight, or obesity), then S can take R = 6 possible values. Each element in S is referred to as a "randomization stratum". Within each randomization stratum, permuted blocks are used for sequential treatment allocation. Each permuted block contains fraction π j j's for j = 0, . . . , J, with j representing treatment group j. At the onset of treatment allocation, a permuted block is randomly chosen and used to sequentially assign treatment.
After a block is exhausted, a new block is used.
Compared with simple randomization, stratified randomization is able to achieve balance within each randomization stratum, i.e., exact fraction π j of participants are assigned to group j. Under stratified randomization, the treatment assignments (A 1 , . . . , A n ) are not independent of each other; and (A 1 , . . . , A n ) are conditionally independent of (W 1 , . . . , W n ) given (S 1 , . . . , S n ). For each participant i, we assume that S i is encoded as R − 1 dummy variables (dropping one level to avoid collinearity) in the baseline vector X i . We finish this section by introducing a few additional definitions. For any two symmetric matrices V 1 and V 2 with the same dimension, we denote V 1 For any estimator ∆ of ∆ * , we call V the asymptotic covariance matrix of ∆ if √ n( ∆ 1 − ∆ * ) weakly converges to a multivariate normal distribution with mean 0 and covariance matrix V. If two estimators ∆ 1 and ∆ 2 of ∆ * have asymptotic covariance matrix V 1 and V 2 respectively, we call ∆ 2 is (asymptotically) equally or more precise than ∆ 1 if Such an expression is commonly used for scalar estimators, and we extend it to vector estimators.

Estimators
For estimating the average treatment effects ∆ * , we first introduce three commonly-used estimators, the ANCOVA estimator (Section 4.1), the MMRM estimator (Section 4.2), and the MMRM estimator with visit-by-covariates interactions (Section 4.2). Next, we propose a new estimator, the IMMRM estimator (Section 4.3), which models treatment heterogeneity and heteroscedasticity based on a variant of the MMRM model. Our main result in Section 5 is that, given Assumption 1 and mild regularity conditions, the IMMRM estimator is equally or more precise than the ANCOVA estimator, MMRM estimator and MMRM-VCI estimator under simple or stratified randomization, without assuming the IMMRM working model is correctly specified.

The Analysis of Covariance (ANCOVA) estimator
The ANCOVA estimator for ∆ * is defined as the maximum likelihood estimator (MLE) for parameters (β A1 , . . . , β AJ ) in the working model below: where (β 0 , β A1 , . . . , β AJ , β X ) are parameters, and ε is independent of X and follows a normal distribution with mean 0 and unknown variance σ 2 . We denote the ANCOVA estimator as ∆ (ANCOVA) .
Since it was first proposed by Fisher et al. (1937), ANCOVA has been extensively studied and applied. Under simple randomization, Yang and Tsiatis (2001) showed that the ANCOVA estimator for two arms (i.e., J = 1 in our setting) is consistent and asymptotically normal given arbitrary misspecification of its working model. Ye et al. (2021)

The mixed-effects model for repeatedly measured outcomes (MMRM)
The MMRM working model is defined as, for each t = 1, . . . , K, where (β 0t , β A1t , . . . , β AJt ) are parameters that are specific for visit t, β X are parameters that are invariant across t, ε = (ε 1 , . . . , ε K ) is independent of X and has a multivariate normal distribution with mean 0 and unknown covariance Σ. The covariance matrix Σ is assumed to be unstructured; that is, no other assumption is made on Σ except that it is positive definite. The MMRM estimator ∆ (MMRM) for ∆ * is defined as the MLE for parameters (β A1K , . . . , β AJK ) in the working model (2) using observed data (O 1 , . . . , O n ).
The MMRM model (2) is derived from a linear mixed-effects model with fixed effects and random effects, where parameters (β 0t , β A1t , . . . , β AJt ) and β X are fixed effects. The random effects are marginalized and implicitly represented in the covariance matrix Σ.
In some clinical trials, e.g. Sorli et al.
which differs from the MMRM working model only on the regression coefficients of X, where β Xt is substituted for β X . We define the MMRM-VCI estimator ∆ (MMRM-VCI) as the MLE for (β A1K , . . . , β AJK ) in model (3). Compared with the MMRM model, the MMRM-VCI model can capture the time-varying correlation of covariates and outcomes. In addition, MMRM-VCI is also a generalization of the constrained longitudinal data analysis by Liang and Zeger (2000), which focused on X being the baseline value of the outcome and following a normal distribution. As we will show in Section 5.3, under two-arm equal randomization, the MMRM-VCI estimator is equally or more precise than the MMRM estimator and the AN-COVA estimator.
In the MMRM model (2) and the MMRM-VCI model (3), the correlation of covariates X and the outcome vector Y is assumed to be not varied among treatment groups, which we call the homogeneity assumption. In addition, the covariance matrix of Y is also assumed to be the same across treatment groups, which we refer to as the homoscedasticity assumption.
Due to these assumptions, although MMRM (or MMRM-VCI) utilizes information from intermediate outcomes, it does not provide guaranteed precision gain compared to ANCOVA if the MMRM (or MMRM-VCI) working model is misspecified. Schuler (2021) showed that, under two-arm equal randomization, ANCOVA is asymptotically more powerful than MMRM when no outcomes are missing. In two of our data applications below (Section 7), the MMRM estimator is less precise than ANCOVA.
Compared with MMRM, the IMMRM working model has two improvements. First, the inclusion of treatment-covariates-visits three-way interaction terms allows the relationship between the outcomes and baseline variables to vary across treatment groups and visits. Such interaction terms models heterogeneity, which has been shown by Tsiatis (2007); Ye et al. (2021) as an effective method to improve precision for scalar outcomes. We extend this idea to longitudinal repeated measures data. Second, the covariance matrix of Y is modeled separately for each treatment group, which accounts for heteroscedasticity. Gosho and Maruo (2018) first proposed the idea of modeling heteroscedasticity in MMRM; however, they only provide empirical results to show its benefits. We show, in Section 5.2, that modeling heteroscedasticity is necessary for achieving asymptotic precision gain when repeated measure outcomes are jointly modeled. We also give an example (in the Supplementary Material) that MMRM-VCI, which does not account for heteroscedasticity, is 5% less precise than ANCOVA, even though MMRM-VCI uses information from intermediate outcomes.
In addition, we provide the conditions for V (est) = V (IMMRM) in the Supplementary Material. For performing hypothesis testing and constructing confidence intervals, we provide consistent estimators for the asymptotic covariance matrices V (est) and V (est) in the Supplemen- (for stratified randomization). In addition, the asymptotic result for the ANCOVA estimator is not new, which has been proved by Yang and Tsiatis (2001); Tsiatis (2007); Ye et al.
(2021); we include it in Theorem 1 for conveniently presenting the comparison of asymptotic covariance matrices. The major innovation, also the major challenge, of the proof is to derive the partial order (5), where the multivariate non-missing indicator M adds substantial algebraic difficulty. We overcome this challenge by developing a series of inequalities related to functions of M and positive definite matrices.

How adjustment for intermediate outcomes improves precision
Consider the IMMRM working model (4) with a modification that intermediate outcomes are excluded, i.e., Then the IMMRM estimator for ∆ * under the above working model (6) is equivalent to the "ANHECOVA" estimator, which is proposed by Ye et al. (2021). As a special case of the IMMRM model (4), the working model (6) ] ⊤ is the covariance between any random vectors U 1 , U 2 with finite second moments, and V ar(X) = Cov(X, X) is the covariance matrix of X.  Unlike the IMMRM estimator, adjusting for the intermediate outcomes as in the MMRM working model (2) or the MMRM-VCI working model (3) may increase the variance, com-pared to the ANCOVA estimator. For the MMRM estimator, the effect of X on Y t (j) is modeled as a constant vector across visits and treatment groups, while the ANCOVA estimator focuses on the last visit and models the effect of X on Y K (j) as a constant vector across treatment groups. If ANCOVA happens to capture the true effect of X on Y t (j), then it is more precise than MMRM, and vice versa. Our data application in Section 7 also shows that MMRM can be either more or less precise than ANCOVA on different data sets.
For the MMRM-VCI estimator, although it models Y K on (A, X) in the same way as in ANCOVA, it fails to capture the heteroscedasticity. In the supplementary material, we give a simple example showing that the MMRM-VCI estimator has a 5% larger asymptotic variance than the ANCOVA estimator in the presence of heteroscedasticity. An exception, which we discuss in detail in Section 5.3, is the two-arm equal randomization, where MMRM-VCI is equally or more precise than ANCOVA.

Special case: two-arm equal randomization
In a general setting, e.g. multi-arm trials or unequal randomization, none of the ANCOVA, MMRM or the MMRM-VCI estimator is asymptotically more precise than the others. Under the two-arm equal randomization (i.e. J = 1 and π 1 = π 0 ), however, the following corollary implies that the MMRM-VCI estimator has equal or smaller asymptotic variance than ANCOVA and MMRM; in addition, the MMRM-VCI estimator has the same asymptotic variance under simple or stratified randomization.
Corollary 2. Assume J = 1, π 1 = π 0 , Assumption 1, and regularity conditions in the Corollary 2 extends the results of Schuler (2021), which assumes no missing data and shows V (MMRM) ≥ V (ANCOVA) = V (MMRM-VCI) under two-arm, simple, and equal randomization. Despite the advantage of MMRM-VCI over ANCOVA and MMRM, the IMMRM estimator remains to be equally or more precise than MMRM-VCI, for which we demonstrate by the data application.
6 Simulation study

Simulation settings
We conducted a simulation study assessing the performance of the ANCOVA (1), MMRM (2), MMRM-VCI (3) and IMMRM (4) estimators in varied simulation settings. This simulation study was based on the IMAGINE-2 study introduced in Section 2.1.
In the simulation, we considered five post-randomization visits (K = 5), three treatment groups (J = 2) and four randomization strata (R = 4), trying to duplicate the setting of the IMAGINE-2 study while considering a multi-arm study. The simulated data were generated by the following steps.
First, we defined the potential outcome for a reference group. We took the 844 completers at the end of week 52 from insulin peglispro arm to serve as the super-population for a reference group, called AC (active comparator). By doing so, the underlying data generating mechanism remained unknown and also mimicked the real data distribution. Let us use Y i (0) to represent the outcome vector for patient i had this patient taken treatment AC.
We then defined the potential outcome for two other treatment groups, named TRT1 and TRT2, as two hypothetical basal insulin treatments. For each i = 1, . . . , 844, we used Y i (1) and Y i (2) to denote the outcome vector for patient i if he/she had been assigned treatment TRT1 and TRT2, respectively. Let X 1i be the baseline HbA1c and X 2i be the baseline indicator of LDL cholesterol ≥ 2.6mmol/L (100mg/dL) for patient i. We generated the potential outcome for Y i (1) and Y i (2) by: where, for j = 1, 2 and t = 1, . . . , K, c t (j) is a constant specifying the average treatment effect comparing TRTa to AC at time t, α t (j) and γ t (j) are coefficients that determine the degree of heterogeneity among treatment arms, and X 1i and X 2 1i are averages of X 1i and X 2 1i across i in the AC group. The quadratic terms adds another layer of model misspecification. We set These values were set to explore a range of mild treatment heterogeneity. The negative signs of α(j) and γ(j) indicated that a higher baseline HbA1c is associated with a larger HbA1c change.
After we defined the potential outcomes for all treatment groups, we simulated the randomized clinical trial by resampling with replacement from the empirical distribution of We considered two settings for sample size: n = 50 or n = 200 per arm. These numbers represent typical phase 2 or 3 diabetes clinical trials. We defined a strata variable S i which covers all the joint levels of stratum defined by baseline HbA1c (I{X i1 ≤ 8.5}) and LDL cholesterol (I{X i2 = 1}). We applied stratified permute block randomization with a block size of 6 to randomly assign the resampled patients to three treatment arms with 1:1:1 randomization ratio. Then the treatment variable was A i taking values in {0, 1, 2} and the outcome vector was Y i = 2 j=0 I{A i = j}Y i (j). In the next step, we generated missing outcomes to the simulated data. We considered two missing data mechanism: missing completely at random (MCAR) and missing at random (MAR), both under the assumption of monotone censoring. We mimicked the missing data percentages across 5 post-baseline visits in IMAGINE-2 study such that 3%, 6%, 10%, 13%, and 15% are expected to be missing at visits 1-5 respectively. For MCAR, the censoring time was generated by a logistic regression (with an intercept only) to achieve the missing data percentages above. For MAR, the censoring time was determined by a logistic regression model on the treatment arms, HbA1c values at the previous visit and baseline. Details of the missing data mechanism under MAR is given in the Supplementary Material. The dropout rate was made higher in the arms AC and TRT1 compared to the arm TRT2 and the patients are more likely to drop out given a higher HbA1c observed from the previous visit.
Given a simulated data, ANCOVA (1), MMRM (2), MMRM-VCI (3) and IMMRM (4) were then used to compute the estimate of the average treatment effect in change from baseline of HbA1c comparing arms TRT1 and TRT2 to arm AC at the last visit, their standard errors and the 95% confidence intervals. For all estimators, the standard error is calculated by the sandwich variance estimator (using the option "empirical" in SAS with an adjustment for the variability in the mean covariates as pointed out by Qu and Luo, 2015 in estimating the average treatment effect). Such estimators do not account for stratified randomization.
The above procedure were repeated for 10,000 times for each simulation setting, i.e.

Simulation results
The simulation results are summarized in Table 1 for n = 50 and Table 2 for n = 200. Across the simulation settings, all candidate estimators provided small bias, coverage probability close to 95% and type-1-error (i.e., the power comparing TRT1 vs AC) close to 5%. Under MAR, all estimators tended to have slightly larger bias than those under MCAR, since none of the estimators are unbiased under MAR and model misspecification. Such bias is related to the magnitude of treatment effect and the missing mechanism.
Among all estimators, the IMMRM estimator had the smallest mean squared error reflected by the RMSE measure. This is consistent with our asymptotic results, and the simulation also demonstrated the advantage of IMMRM when the sample size is small. The MMRM estimator, however, had the largest mean squared error among the four estimators.
The MMRM-VCI estimator showed a precision gain compared to ANCOVA and MMRM, while such gain is not asymptotically guaranteed as discussed in Section 5.2.
Compared with empirical standard error and RMSE, the average standard error, which reflected the current practice, failed to account for the precision gain from the stratified randomization. For ANCOVA, MMRM and MMRM-VCI, the average standard error were in general greater than the empirical standard error, while the IMMRM estimator appeared to have similar average and empirical standard errors, which further supported Theorem 1.
For MMRM estimator, ignoring stratified randomization led to an overestimation in standard error by up to 25%.
All three trials used stratified randomization, which we accounted for in the computation of standard error by using the consistent variance estimator given in the Supplementary Material. For each estimator and each treatment comparison, we computed the estimate, the standard error (SE) and the proportional variance reduction (PVR), defined as one minus the variance ratio of the estimator and the MMRM estimator. Positive PVR is in the direction of variance reduction, while negative PVR indicates that the compared estimator has larger variance than the MMRM estimator. The PVR is a measurement of precision change of using the estimator compared with the MMRM estimator, and can be translated into the same amount of proportional sample size change for achieving a desired power.
For Trial 1, due to very limited data for participants with no sulfonylurea/meglitinide use at baseline, we drop this stratification variable in the analysis, resulting in 4 randomization strata in total. Combining small strata, as we did here, is a method that achieves better finite sample estimation and controls the type I error at the same time, as discussed by Wang et al. (2019b). Similar to Trial 1, we dropped the stratification variable sulfonylurea/meglitinide use from Trial 2 due to its data limitation, resulting in 12 total strata. For Trial 3, we dropped the metformin use variable due to limited data for participants with no metformin use, resulting in 4 randomization strata in total. Table 3 summarizes our data applications, which consists of six treatment-control comparisons. Among all six comparisons, the IMMRM estimator had the smallest standard error in five. All the estimators provided similar estimates. In the comparison (tzp 1mg vs pbo) that IMMRM was less precise than MMRM-VCI, it still outperformed ANCOVA and MMRM; such results may occur in practice when the sample size is small, the intermediate outcomes are prognostic, and the homogeneity and homoscedasticity assumptions hold.
Overall, the IMMRM estimator was 2-24% more precise than ANCOVA, 5-16% more precise than MMRM, and up to 11% better than the MMRM-VCI. Besides IMMRM, MMRM-VCI also showed precision gain compared to standard practice that uses ANCOVA or MMRM.
Although such gain is not guaranteed, MMRM-VCI may have comparable standard error compared to IMMRM in practice. Table 3 also shows that MMRM may have a lower or higher variance than ANCOVA, which indicates that adjustment for intermediate outcomes using MMRM may harm the precision. In Trial 2, the precision loss was 14% compared to ANCOVA. We hence recommend using IMMRM or MMRM-VCI, instead of MMRM, to adjust for intermediate outcomes.
While Trials 1, 2 and 3 are different in many aspects, the IMMRM estimator had smaller standard error than the ANCOVA and MMRM estimators in all three trials. Trial 1 had a large sample size, unequal two-arm randomization and five post-randomization outcomes.
Trial 2 used two-arm equal randomization with three post-randomization outcomes. In addition, trial 2 had small percentage of missing data in the primary outcome (less than 6%) across all visits. Our analyses for Trial 3 involved five treatment groups, eight postrandomization outcomes and 19% missing outcomes at the last visit, which is a typical setting of phase 2 trials. By the data application on the three clinical trials, we were able to evaluate the performance of estimators in three distinct live scenarios with different outcome and missing data distributions, phases of the study, sample sizes, and numbers of treatment groups and clinical visits.

Discussion
For the analysis of longitudinal repeated measures data, we propose the IMMRM estimator that can improve precision, compared to standard practice that uses ANCOVA or MMRM.
In some cases, the variance reduction can be substantial, e.g., 16% in our data application of Trial 2 compared to MMRM. Such precision gain comes from modeling treatment heterogeneity and appropriately adjusting for intermediate outcomes, and can be translated into sample size reduction when planning a trial.
We recommend applying the IMMRM estimator for large trials, which can lead to more sample size reduction (Wang et al., 2019a). In addition, since the IMMRM working model involves substantially more parameters than MMRM (K(p+1)(J +1) parameters in IMMRM compared with p + K(J + 1) parameters in MMRM), large trials, compared to small trials, have more degrees of freedom for the variance estimation. Under the two-arm equal randomization, the MMRM-VCI estimator is a good alternative to the IMMRM estimator for small trials, since MMRM-VCI involves fewer parameters and still outperforms ANCOVA and MMRM.
Leveraging post-randomization information to improve precision is not limited to intermediate outcomes. Our main results also apply to adjustment for other post-randomization continuous-valued random variables measured before the final outcomes, such as the body mass index measured at each visit. When these additional variables provide new prognos-tic information beyond intermediate outcomes and covariates, adding them to the IMMRM model can bring further asymptotic precision gain. A trade-off in finite sample variance estimation, however, is the improved precision versus the reduced degrees of freedom for variance estimation due to estimating additional parameters.
We assumed missing completely at random, which is a typical assumption for proving model-robustness. In contrast, under the missing at random assumption (MAR), i.e. M t (j) Among all variants of MMRM discussed in this paper, the random effects are assumed to be independent of covariates. If, otherwise, random effects are modeled to involve covariates, our main results may not apply. Under simple randomization and correctly specified MMRM model, Cnaan et al. (1997) stated that the MMRM estimator is consistent and asymptotically normal, whose results can be directly generalized to the MMRM-VCI and    Supplementary Material to Robustly leveraging the post-randomization information to improve precision in the analyses of randomized clinical trials In Section 1, we give notations used in later sections. In Section 2, we provide the regularity conditions for Theorem 1 of the main paper. In Section 3, we define the variance estimators for V (est) and V (est) . In Section 4, we introduce a few lemmas for proving our main results. In Section 5, we prove Theorem 1, Corollary 1 and Corollary 2 of the main paper. In Section 6 6, we give an example where MMRM-VCI is less precise than ANCOVA.
In Section 7, we provide the missing data mechanism for MAR in the simulation study.

Notations
Let 1 K be the column vector of length K with each component equal to 1, 0 K be the column vector of length K with each component equal to 0, and e t be the column vector of length K with the t-th entry 1 and the rest 0. Let I K be the K × K identity matrix. Let ⊗ be the , 1}, t = 1, . . . K}. For any matrix (or vector), we use || · || to denote its L 2 matrix (or vector) norm. For any vector v = (v 1 , . . . , v L ), we use diag{v} or diag{v l : l = 1, . . . , L} to denote an L × L diagonal matrix with the diagonal entries being (v 1 , . . . , v L ). For any sequence x 1 , . . . , x n , . . . , we define P n x = n −1 n i=1 x i .

Variance estimators in Theorem 1
For an M-estimator θ ∈ R q defined by P n ψ(A, X, Y , M ; θ) = 0, its sandwich variance estimator under simple randomization is defined as Since ∆ is embedded in θ, we can find C ∈ R J×q such that ∆ = Cθ. Then the variance estimator of ∆ under simple randomization is defined as C V n (ψ, θ)C ⊤ .
We next define V (est) n for est = ANCOVA, MMRM, MMRM-VCI. Define where β X is the MLE of β X in the MMRM working model (2).
Following Equation (16), we define Following Equation (15), we define Similarly, we define  (3) Let A ∈ R K×K be a positive definite matrix. Then

The equality holds if and only if
(4) Let B ∈ R K×K be a positive semi-definite matrix. Then

The equality holds if and only if
Proof. (1) which implies that E[V M ] is positive definite and so invertible.
(2) Consider the case K = 1. we have that Σ reduces to a positive number σ and V m (σ) = where D m ∈ R (K−1)×( K t=1 mt−1) is a matrix taking the first K−1 rows and first ( K Using the formula of block matrix inversion, we can compute that, if m ∈ Ω K , then Since we have shown that E[V M (Σ)] is positive definite and P (M K = 1) = p e K + m∈Ω K p m , by using the formula of block matrix inversion again, we have Since G −1 Ω K is positive semi-definite, then the above derivation shows that the equality holds if and only if where σ j,K = e ⊤ j Σe K , which completes the proof.
. We have the following derivation: where the last equation comes from the fact that The equality holds if and only if V ar(x ⊤ M e t ) = 0 for t = 1, . . . , K, which is equivalent to x ⊤ M being a constant vector.
(5) Define B = C − Σ. The statement is proved by the following derivation: where the first inequality results from Lemma 1 (3), and the second inequality comes from Lemma 1 (4).
For each i, let n i = K t=1 M it be the number of non-missing outcomes and Y o i ∈ R n i be the observed outcomes if n i > 0. Let t i,1 < · · · < t i,n i denote the ordered list of visits when the outcomes are not missing. For example, t i,1 is the first non-missing visit for subject i.

The corresponding log likelihood function conditional on
which implies that ∆ (MMRM) and ∆ (MMRM-VCI) are M-estimators. In the above expression of ψ (MMRM) , we omit (Σ) from V M (Σ) for conciseness. We note that V M is a random matrix taking values in R K×K and defined in the same way as in Lemma 1.
Following a similar procedure as for the MMRM model (2), the estimating functions for the MLE θ (IMMRM) under the MMRM model (3) are where Lemma 4,and IF (est) represents the J-dimensional influence function. We note that Theorem 1 of Wang et al.
(2019) is developed for binary treatment (i.e. J = 1) and scalar outcome (i.e. K = 1), but their proof can be easily generalized to accommodate multiple treatment arms and repeated measured outcomes (as in Example 3 of Wang et al., 2019).
Following a similar procedure, we next show that ∆ (IMMRM) is model-robust. We have that, for j = 1, . . . , J, being the probability limit of Σ j ( α j ) in the IMMRM model (3) for j = 0, . . . , J. The assumption P (M (j) = 1 K ) > 0 and Lemma 1 (1) implies that E[V jM ] is invertible. Thus, for We next prove that √ n( ∆ (est) − ∆ * ) weakly converges to a normal distribution, under simple or stratified randomization. Given Equations (6), it suffices to show that . . , n are independent to each other and identically distributed. Since the regularity conditions implied that IF (est) has finite second moment, then the central limit theorem implies the desired weak convergence. Furthermore, we have By the regularity conditions, we apply Lemma 3 and get which completes the proof of asymptotic normality. In addition, Lemma 3 also implies that For the ANCOVA, MMRM, MMRM-VCI and IMMRM estimators, the influence func-tions by Lemmas 5 and 6 are given below: where Furthermore, we have where We next compute V (est) . Lemma 6 implies that, for j = 0, . . . , J, ). Then using Equation (8), we get where e j+1 ∈ R J+1 has the (j + 1)-th entry 1 and the rest 0. Hence by Lemma 3, we have where where is positive semi-definite, and Using Equations (11) and (14) where π = (π 0 , . . . , π J ) ⊤ and In the above derivation, the first " " results from Lemma 1 (2), the second " " comes from P * (M K = 1) ≤ 1 and Lz ⊤ = Lr ⊤ . Since diag{π} ππ ⊤ , then (diag{π} − ππ ⊤ ) ⊗ I p is positive semi-definite (by Lemma 2) and hence V (ANCOVA) V (IMMRM) .
Next, for showing V (MMRM-VCI) V (IMMRM) , we can follow a similar proof as in the previous paragraph, where Σ (MMRM) is substituted by Σ (MMRM-VCI) , and get which is positive semi-definite.

Condition (ii) is equivalent to
we get that the (j, j)-th entry of matrix LU S=s (diag{π} − ππ ⊤ )U S=s L ⊤ is which is equal to 0 if and only if either u 0 = u j = 0, or π 0 + π j = 1 and π 0 u j + π j u 0 = 0.
Since we have shown the model-robustness of ∆ (IMMRM) in the proof of Theorem 1, then which implies U 2 = Lr ⊤ (X − E[X]).

Proof of Corollary 1
Proof. The ANHECOVA estimator is a special case of the IMMRM model setting K = 1.