Design and analysis of individually randomized multiple baseline factorial trials

TL;DR

This paper introduces a new trial design for evaluating multiple behavioral interventions in rare diseases with small sample sizes.

Contribution

The novel individually randomized multiple baseline factorial design (MBFD) allows efficient evaluation of multiple interventions with limited participants.

Findings

MBFD requires fewer participants while maintaining statistical power for evaluating multiple interventions.

GEE is recommended over LMM to handle random effect misspecifications in small samples.

Small sample corrections like Mancl and DeRouen estimator improve accuracy for sample sizes below 120.

Abstract

Assessing the effectiveness of behavioral interventions in rare diseases is challenging due to extremely limited sample sizes and ethical challenges with withholding intervention when limited treatment options are available. The multiple baseline design (MBD) is commonly used in behavioral science to assess interventions, while allowing all individuals to receive the intervention. MBD is primarily used to evaluate a single intervention so an alternative strategy is needed when evaluating more than one intervention. In this case, a factorial design may be recommended, but a standard factorial design may not be feasible in rare diseases due to extremely limited sample sizes. To address this challenge, we propose the individually randomized multiple baseline factorial design (MBFD), which requires fewer participants but can attain sufficient statistical power for evaluating at least two interventions and their combinations. Furthermore, by incorporating randomization, we enhance the internal validity of the design. This study describes the design characteristics of a standard MBFD, clarifies estimands, and introduces three statistical models under different assumptions. Through simulations, we analyze data from MBFD using linear mixed effect models (LMM) and generalized estimating equations (GEE) to compare biases, sizes, and power of detecting the main effects from the models. We recommend using GEE to mitigate potential random effect misspecifications and suggest small sample corrections, such as Mancl and DeRouen variance estimator, for sample sizes below 120.

Figures

A standard multiple baseline factorial design in which individuals start in the control condition (*light green*) and then transition into either intervention A (*yellow*) or B (*blue*), depending on their randomization. Individuals spend one time interval receiving the single intervention, then transition into the combined intervention (*dark green*). Measurements are taken repeatedly at each time interval

Full text

Introduction

Conducting clinical trials for rare diseases presents unique challenges that complicate the traditional trial design and statistical analysis. One major difficulty is the small sample size, as rare diseases inherently affect a limited population, making it difficult to recruit enough participants to achieve sufficient statistical power (Partington et al., 2022). Additionally, in some cases where low-risk interventions are being tested, researchers may prefer a trial design in which all participants receive every intervention, rather than being assigned to separate treatment arms, to maximize the available data, account for individual variability, and allow access to potentially effective interventions in areas of high therapeutic need. These constraints necessitate innovative trial designs to ensure robust and reliable conclusions.

To address these challenges, researchers can turn to methodologies developed for rigorous inference in small populations, such as single-case experimental designs (SCEDs). SCEDs represent a family of designs that use repeated measurements to establish a stable baseline and demonstrate experimental control after an intervention is introduced, allowing for valid causal inferences even with a single participant (Epstein & Dallery, 2022; Krasny-Pacini & Evans, 2018; Smith, 2012). The multiple baseline design (MBD) is one of the most frequently used SCEDs in social and behavioral sciences (Baer et al., 1968; Coon & Rapp, 2018). In a standard MBD, all individuals begin in a control or SoC condition, and the intervention is then introduced at staggered times across different sequences (Watson & Workman, 1981). The timing of receiving intervention can either be randomized (experimental design) or non-randomized (quasi-experimental design) (Kratochwill & Levin, 2010; Levin & Ferron, 2021). However, when it is possible, randomization in MBD enhances the internal validity, reduces the potential confounding, and ensures that characteristics across the randomization units are more comparable across different baseline sequences.

Usually, individuals in an MBD are followed up over time, and multiple outcome measures are taken before and after the intervention to establish a stable estimate of the intervention effects. This approach is an extension of the simple before-and-after design and is a robust way of assessing the effectiveness of the intervention as it rules out the potential impact of concurrent events that may be confounded with the intervention (Epstein & Dallery, 2022; Hawkins et al., 2007). This design involves both within- and between-sequence comparisons (Kennedy, 2022). An MBD assumes that if similar effects are observed before and after treatment initiation in different sequences, then this effect of treatment is unlikely to be a coincidence (Slocum et al., 2022).

When randomization is used, MBD may have several advantages over the conventional parallel-arm trial design, which compares an intervention against the SoC. Firstly, individuals receive both the SoC and the intervention, critical to avoid ethical issues around withholding a low-risk, potentially beneficial intervention (Binik, 2019). An MBD design may also boost recruitment as individuals who are interested in the intervention may only participate if they are assured of receiving the intervention at some point. Critically for rare diseases, MBD may also be the only design that provides sufficient statistical power given the available sample size (Sundin & Crespi, 2022).

Sometimes, investigators are interested in testing multiple interventions simultaneously against a SoC (Dziak et al., 2012). In this case, a factorial design is a common choice (Cipriani & Barbui, 2013). However, a standard factorial design does not share the same advantages as an MBD and requires a strong assumption that the interventions are independent (i.e., no interaction effects) to achieve maximum efficiency. Therefore, a standard factorial design may not always be ideal or even feasible, especially when the underlying disease has low prevalence.

A design that combines features of both factorial and MBD is available for evaluating cluster-level interventions, known as a stepped wedge factorial design (Hemming et al., 2015; Hussey & Hughes, 2007; Lyons et al., 2017; Ouyang et al., 2022). However, a factorial MBD that includes individual-level randomization and interventions has not been presented. Therefore, we propose a novel individually randomized multiple baseline factorial design (MBFD) to extend these ideas for individual-level interventions.

As a motivating example, deep vein thrombosis (DVT) is rare but the most common type of venous thrombosis in children (O’Brien et al., 2022; Raffini et al., 2009). The most common chronic complication of DVT is post-thrombotic syndrome (PTS) (M. L. Avila et al., 2016), which is a form of chronic venous insufficiency that manifests with edema, pain, and poor endurance in the affected limb. Compression therapy, considered as SoC, is the only management strategy for pediatric PTS (Mutlak et al., 2019). While it is effective at reducing PTS severity (L. Avila et al., 2021), poor adherence has been a concern (Avila et al., 2021; Montoya et al., 2016). Therefore, a strategy to improve adherence is needed. Compression therapy may also be unsuitable for some children, with one-third of the children finding compression garments too uncomfortable or too difficult to put on, or that their symptoms worsened with compression (Avila et al., 2021).

Drug interventions are gaining popularity for adults but are currently not used in pediatrics (Avila et al., 2024). Thus, trials of these medical interventions are needed in children. Since both an adherence-enhanced compression garment wear program and the use of pharmacological management may benefit children, we want to design a randomized trial to investigate the efficacy of individual and combined programs to decrease PTS symptom severity, measured by the CAPTSure score^©^ (Avila et al., 2019), a validated outcome measure in pediatric patients with lower limb PTS.

Since pediatric DVT is rare, there are inevitable recruitment challenges with a maximum recruitment for the study of around 30 patients. Given this restriction, the standard factorial trial would not provide sufficient statistical power to detect a treatment effect. Participants may also be unwilling to participate unless they are assured of receiving both new interventions at some point. Thus, we developed an individually randomized MBFD to address these two challenges and provide reasonable statistical power.

We critically examine and extend the design and analysis framework of MBFDs, with particular attention to their application in limited sample size settings such as rare diseases or early-phase interventions. We begin by formally defining the MBFD structure, emphasizing its deviation from traditional factorial and multiple baseline designs. We then clarify the estimands that can be unbiasedly identified under its constraints. Given the sequential exposure to treatments and the potential confounding from treatment order and carryover, we carefully consider the causal assumptions required to interpret effects, particularly when estimating combined treatment effects. We explore several statistical models tailored to this design, outlining their assumptions, limitations, and alignment with the estimands. A simulation study evaluates model performance under various conditions, including realistic intraclass correlation coefficients and empirically motivated treatment effect sizes, addressing concerns about generalizability and robustness. Finally, we discuss its broader implications, limitations, and avenues for future methodological development, particularly in improving causal interpretability and estimation efficiency.

Design characteristics

Standard multiple baseline factorial design

Figure 1 introduces the standard MBFD to investigate the impact of two interventions (A and B) and their combination, which can be extended to more than two interventions. In this design, the outcome is repeatedly measured, with the number of outcome measurements and the time between outcome measurements selected to be sufficient to understand the variability in the outcome. Each column in Fig. 1 represents a time period. Each box in Fig. 1 represents a time interval in which the outcome is measured. The MBFD includes multiple treatment sequences that specify the number of intervals for which the patients will receive the control condition, the initial intervention, and the combination. We recommend a minimum of three sequences per intervention (Lanovaz & Turgeon, 2020), with everyone randomized to one of the sequences (rows in Fig. 1). After randomization, everyone begins the study receiving the control condition (light green cells). Following the initial treatment phase, individuals transition into either intervention A or B, and remain with the selected intervention for one interval before transitioning into the combined intervention. There are at least five intervals, and measurements are taken at the end of each interval for every individual. The standard MBFD design also ensures that everyone receives the combined intervention for one interval. Fig. 1A standard multiple baseline factorial design in which individuals start in the control condition (light green) and then transition into either intervention A (yellow) or B (blue), depending on their randomization. Individuals spend one time interval receiving the single intervention, then transition into the combined intervention (dark green). Measurements are taken repeatedly at each time interval

In our proposed MBFD, the combined intervention (AB) is always delivered sequentially following exposure to either A or B. As such, the AB effect should be interpreted as a conditional sequence effect. That is, the impact of receiving both interventions in sequence compared to SoC. This differs from the concurrent AB effect estimated in a traditional factorial trial. This sequential approach is particularly relevant in rare-disease settings where the sample size required for a pure factorial design is unattainable, and where maximizing exposure to potentially beneficial treatments is desirable. Our design differs from a conventional four-arm factorial design, in which four concurrent groups (e.g., control/SoC, A, B, and AB) are compared, and AB represents simultaneous delivery of both interventions. Such a structure allows direct estimation of main effects and the unconditional AB interaction. This choice to include a sequential combination was guided by the limited sample size available for rare disease trials and the desire to preserve the strengths of the MBD framework while exploring the effects of combined interventions.

Due to these limitations, the objective of our motivating trial is, therefore, to test the following three hypotheses:

1. A behavioral intervention to improve adherence (A) alone decreases PTS symptom severity compared to SoC.
2. Pharmacologic management of PTS (B) decreases symptom severity compared to SoC.
3. The sequential combination of A and B, following prior exposure to one of them, decreases PTS symptom severity compared to SoC.

Extensions to standard multiple baseline factorial design

In Fig. 1, we demonstrate a standard design where all individuals receive the single intervention for a single time interval before transitioning to the combined intervention. This design can be modified for different applications. For example, if more than one outcome measurement is required to show intervention effects, then the number of intervals for the intervention can be extended. Secondly, the standard design assumes that all individuals exit the trial at the same time, which is convenient but individuals receive the combined intervention for a varied number of intervals. This can also be modified so that individuals leave the trial after receiving the combined intervention for a sufficiently long time. The number of sequences chosen is a crucial aspect of MBDs. The inclusion of more sequences can enhance the reliability of the conclusions (Lanovaz & Turgeon, 2020) and studies have also suggested increasing the length of time between the start of the intervention in adjacent sequences. For example, instead of having the intervention rolled out every time interval, this can be increased to three or more intervals, which can further reduce the impact of coincidental events, and ensure the observed effects is attributed to the interventions (Slocum et al., 2022).

Statistical model

In this section, we introduce statistical models that can analyze data from an MBFD.

MBFD without interaction effect (Model 1)

In a standard MBFD, each participant is followed over time and their outcomes are measured repeatedly. Let us assume that we have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} participants, and each one contributes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J$$\end{document} observations. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Y}_{ij}$$\end{document} denote the continuous outcome of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i-th$$\end{document} individual at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j$$\end{document} . Under the assumption of conditional independence (given random effects), the data can be modeled as, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Y}_{ij}$$\end{document} can be modeled as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Y}_{ij}=u+ {\beta }_{T}*j+{\beta }_{A}{X}_{Aij}+{\beta }_{B}{X}_{Bij}+{\alpha }_{i}+{\varepsilon }_{ij}$$\end{document}

Here, the overall baseline mean is denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u$$\end{document} . \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{T}$$\end{document} , a continuous time interval effect, is used to control for the underlying time trend. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X}_{Aij}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X}_{Bij}$$\end{document} are binary indicators for whether individual \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i$$\end{document} received intervention A and/or B at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j$$\end{document} , taking the value of one when the individual receives the corresponding intervention; otherwise, it is zero. The random intercept term, which links observations for each individual, is expressed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{i}$$\end{document} , following a normal distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N(0, {\tau }^{2})$$\end{document} . Finally, the residual error term is expressed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varepsilon }_{ij} \sim N(0, {\sigma }^{2})$$\end{document} .

Model 1 assumes no A and B interaction and no history/order (sequence) dependence (i.e., the effect of A does not depend on prior or concurrent exposure to B, and vice-versa). In this longitudinal MBFD, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} are parameters for the averaged treatment effects of interventions A and B over the exposed intervals. MBFD staggers the timing of intervention implementation across different sequences, so the baseline and intervention phases are not contemporaneous across each of the sequences (e.g., the baseline and intervention data are not collected simultaneously). This requires us to control for the potential underlying time trends.

During the combined (sequential) phase, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X}_{Aij}=1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X}_{Bij}=1$$\end{document} ; under the additivity and no sequence and carryover effect assumption, the expected effect of receiving both interventions equals \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} . This assumption may be violated in practice, especially under sequential delivery, so it should be evaluated contextually. If interaction or historical effects are present, estimates of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} , and the implied AB effect will be biased under Model 1.

MBFD with interaction effect (Model 2)

To relax the assumption of additive effects, we introduce an interaction term between the two interventions:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Y}_{ij}=u+ {\beta }_{T}*J+{\beta }_{A}{X}_{Aij}+{\beta }_{B}{X}_{Bij}+{\beta }_{I}{X}_{Aij}{X}_{Bij}+{\alpha }_{i}+{\varepsilon }_{ij}$$\end{document}

In this model, we add a single interaction term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X}_{Aij}{X}_{Bij}$$\end{document} . Because AB is delivered sequentially, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document} quantifies the departure from additivity under sequential delivery (no sequence effects). That is how the effect of the second intervention changes given prior exposure to the first, which should not be interpreted as a concurrent factorial interaction. If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}\ne 0$$\end{document} , the combined treatment effect deviates from the additive sum ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} ). Thus, this model allows estimation of whether the effect of one intervention is modified by the presence of the other in sequential delivery.

Alternative multiple-arm model (Model 3)

As an alternative, the combined intervention can be modeled as a distinct treatment arm:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Y}_{ij}=u+ {\beta }_{T}*J+{\beta }_{A}{Z}_{Aij}+{\beta }_{B}{Z}_{Bij}+{\beta }_{C}{Z}_{Cij}+{\alpha }_{i}+{\varepsilon }_{ij}$$\end{document}

In this specification, binary indicators \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Z}_{Aij}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Z}_{Bij}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Z}_{Cij}$$\end{document} are mutually exclusive binary indicators for A, B, and the combination AB (denoted as C). The coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ​ represents the average effect of sequential exposure to both interventions vs. SoC under the randomization. The incremental effect of adding the second intervention can be obtained by comparing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ​ with the appropriate main effect ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} ​ or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} ​), averaged across randomization order.

Estimands

Specifying estimands is essential to interpreting model parameters and aligning them with the study’s clinical or scientific objectives (Kahan et al., 2024). In MBFDs, sequential exposure complicates the interpretation of treatment effects. Drawing from the factorial design literature (Kahan et al., 2022), we highlight three relevant estimands:

• E1: Effect of intervention A (or B) in the absence of the other.
• E2: Effect of intervention A (or B) in the presence of the other.
• E3: Combined effect of interventions A and B when delivered sequentially, following prior exposure to A or B.

To aid interpretation, we provide a structured comparison of the three statistical models considered in this manuscript. Table 1 summarizes the key assumptions, parameters of interest, and how the combined effect of interventions A and B is represented under each approach. In Model 1, the combined effect is assumed to be additive (​ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} ) and unbiased estimation relies on the absence of interaction or sequence effects. Model 2 relaxes this assumption by including an interaction term, allowing the combined effect to depart from additivity ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document} ). Finally, Model 3 treats the combined intervention as a distinct treatment condition, estimating its effect directly ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ) without requiring additivity. This structured comparison highlights how each model targets a slightly different estimand and underscores the importance of selecting a model consistent with the trial’s objectives and underlying scientific assumptions. Table 2 shows how each of the estimand match our hypotheses for example trial. Table 1. Estimands in MBFD and corresponding estimates from all modelsModelKey AssumptionsParameters of InterestRepresentation of Combined EffectModel 1: Additive (no interaction)- No A×B interaction- No sequence effects- No carryover effect β A: Average effect of A vs SoC

β B: Average effect of B vs SoC

β A + βB Model 2: With interaction- Interaction allowed- No sequence effects β A: Effect of A vs SoC (without B)

β B: Effect of B vs SoC (without A)

β I: Departure from additivity under sequential delivery

β A + βB + βI Model 3: Multi-arm specification- AB treated as a distinct condition- No additivity assumption required- No sequence effects β A: Effect of A-only vs SoC βB: Effect of B-only vs SoC βC: Effect of AB vs SoC

β C(can be interested as the averaged effect across orders if sequence effects exist)Table 2. Estimands in MBFD and corresponding estimates for the hypotheses in example trialsHypothesisCorresponding estimandModel(s) for unbiased estimationKey assumptionsH1: A behavioral intervention to improve adherence alone decreases PTS symptom severity compared to SoCEffect of A alone vs. SoC (main effect of A in absence of B)- Model 1 (Additive): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} estimates effect of A (valid if no interaction or sequence dependence)- Model 2 (Interaction): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} gives effect of A when B = 0 (directly unbiased even if interaction exists)- Model 3 (Multi-arm): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} (from A vs. SoC) directly identifies A’s effect- No carryover beyond modeled terms- For Model 1: assumes no interaction and no sequence effects- For Model 2: unbiased if order effects do not matter (assume on sequence effect)H2: Pharmacologic management of PTS decreases symptom severity compared to SoCEffect of B alone vs. SoC (main effect of B in absence of A)- Model 1 (Additive): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} estimates effect of B (valid if no interaction or sequence dependence)- Model 2 (Interaction): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} gives effect of B when A = 0 (robust to interaction)- Model 3 (Multi-arm): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} (from B vs. SoC) directly estimates B’s effect- For Model 1: assumes independence (no A–B interaction)- For Model 2: unbiased if no sequence effects H3: The sequential combination of A and B, following prior exposure to one of them, decreases PTS symptom severity compared to SoCEffect of A then B (or B then A) vs. SoC (averaged sequential effects)- Model 1 (Additive): E3 = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}+{\beta }_{B}$$\end{document} (valid only if no interaction/sequence effects)- Model 2 (Interaction): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}+{\beta }_{B}+{\beta }_{I}$$\end{document} gives combined sequential effect (robust if order doesn’t matter)- Model 3 (Multi-arm): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} directly estimates the AB sequence effect vs. SoC- For Model 1: requires strict additivity and no sequence dependence- For Model 2: assumes no sequence effects; unbiased if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document} captures modification correctly

Estimation methods

To obtain estimates from an MBFD, we need to model the correlations between the observations from the same individual but different intervals. Usually, data of this type are modeled using either linear mixed effect model (LMM) or generalized estimating equations (GEE) (Dahmen & Ziegler, 2004; Molenberghs & Verbeke, 2000). Both models can be easily fitted using standard statistical software, such as R, SAS, and Stata. Previous studies have shown that, for continuous outcomes, LMM and GEE produce very similar results (Ouyang et al., 2024).

Models 1 to 3 are examples of LMM. Compared to LMM, GEE uses robust variance estimators that provide reliable estimates even when the random effects are not correctly specified (Gardiner et al., 2009; Liang & Zeger, 1986). For instance, Models 1 to 3 assume that there is only a random intercept with an exchangeable correlation structure, which implies that any two observations from the same individual have the same correlation regardless of when they were measured. However, in real applications, this correlation structure may be more complex, such as auto-regressive (AR) or even unstructured. The actual random effects may also have more components. For example, instead of having random intercept only, the actual random effects may contain both random intercept and slope. In this case, LMM could lead to invalid statistical inferences (Ouyang et al., 2023a, 2024), but GEE would still provide valid estimates in general (Ford & Westgate, 2020). However, in GEE, the classic robust variance estimator tends to underestimate the standard error when the sample size is small. Therefore, a small sample correction is often required in these cases (Ford & Westgate, 2020; Westgate & Burchett, 2016).

Simulation study

Estimators for the same quantities from different models may have different statistical properties. Therefore, it is critical to evaluate these properties to make recommendations about the appropriate model to use for MBFD. We conducted a simulation to compare the unbiasedness of estimators, the size, and power of hypothesis tests in the context of the MBFD. Specifically, we compared estimates from the three proposed models (Model 1, 2, and 3) using two different analytical approaches (LMM and GEE).

Data generation

In our simulation study, we modeled data using the standard MBFD structure illustrated in Fig. 1. Participants were individually randomized to one of six intervention sequences. Each individual was observed over five intervals, beginning with a control condition, followed by either intervention A or B, and finally the combined intervention (A + B). The outcome data were generated from Model 3 by assuming that the baseline means ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u$$\end{document} ) was zero, and the linear time effect ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{T}$$\end{document} ) was one. Without loss of generalizability, we assumed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} equals one. To reflect a range of plausible within-individual correlations, we varied the random intercept standard deviation ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} ) to achieve intraclass correlation coefficients (ICC) of 0.05, 0.10, and 0.30. These correspond to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\in \{0.23, 0.33, 0.65\}$$\end{document} , respectively. Three treatment effect scenarios were considered: (1) null effects \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}]={[0.0, 0.0, 0.0]}^{T}$$\end{document} , (2) additive effects (no interaction): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}]={[0.8, 0.8, 1.6]}^{T}$$\end{document} or (3) non-additive effects (with interaction): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${[0.8, 0.8, 2.0]}^{T}$$\end{document} . The magnitude of the effects was deliberately set to be large. This ensures high statistical power in some conditions, allowing for comparisons of power between the LMM, GEE, and GEE-MD models. If the effects were too small, power would be universally low, making it impossible to discern differences in the relative efficiency of the methods.

To examine how model performance varies with study size, we considered three sample sizes: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30, 60, and 120. These values were selected to reflect standard, large, and extremely large sample sizes in rare disease trials. Specifically, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30 represents a sample size for single-site rare disease trials, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 60 reflects a multi-site study for individually randomized MBFDs in rare disease, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 120 provides insight into performance as statistical power increases and estimation becomes more stable. By simulating across this range, we aimed to assess how bias, type I error, and power behave under realistic constraints, as well as under more favorable conditions where larger sample sizes are available. This helps clarify the strengths and limitations of MBFDs and the corresponding statistical models across a spectrum of practical applications.

Estimands

We were primarily interested in the effect of interventions A and B in the absence of the other intervention, which can be directly estimated by the beta coefficients ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} ) from three models. We presented results for the other parameters ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ) for completeness.

Data analysis

We used both LMM and GEE to analyze the simulated trial data. Three LMM, shown as Models 1, 2, and 3, were fitted using “lmerTest” package in R (Kuznetsova et al., 2020). Satterthwaite degree of freedom is used in LMM by default. Additionally, three GEE models were also fitted to each simulated dataset using “geepack” package in R (Højsgaard et al., 2024). Degrees of freedom equal to the number of individuals minus the number of parameters in the model were used. These GEE models have the same marginal component as Models 1, 2, and 3 with an exchangeable correlation structure. Since our sample size is as low as 30, we fitted both standard GEE and GEE with Mancl and DeRouen variance estimator (GEE-MD) to correct the biased estimated standard error in GEE when the sample size is small (Mancl & DeRouen, 2001).

Performance measures

To evaluate the models, we focus on three performance measures, bias, type I error and power.

To estimate bias, we extract the estimated value for parameters of interest and compare it with the true value. Assuming \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} is the parameter of interest, the bias can be estimated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{{n}_{sim}}{\sum }_{i=1}^{{n}_{sim}}(\widehat{{\beta }_{i}}-\beta )$$\end{document} (where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${n}_{sim}$$\end{document} is the number of simulations) (Morris et al., 2019).

To determine the type I error rates, we used the hypothetical trial data simulated under null ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}]={[0, 0, 0]}^{T}$$\end{document} ). This enabled us to test three null hypotheses under three models (Model 1: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\eta}}}_{1}=[{\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}]$$\end{document} = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${[0, 0]}^{T}$$\end{document} ; Model 2: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\eta}}}_{2}=[{\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}]$$\end{document} = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${[0, 0, 0]}^{T}$$\end{document} ; Model 3: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\eta}}}_{3}= [{\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}]$$\end{document} = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${[0, 0, 0]}^{T}$$\end{document} ). Then, we calculated the percentage of simulations that falsely rejected the null hypothesis. We used a type I error rate of 0.05 and Bonferroni corrections were applied (e.g., 0.05/2 for Model 1 and 0.05/3 for Models 2 and 3).

To determine power, we used hypothetical trial data simulated under alternative hypothesis,either \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{\beta }_{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}]={[0.8, 0.8, 1.6]}^{T}$$\end{document} (without interaction effect) or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${[0.8, 0.8, 2.0]}^{T}$$\end{document} (with interaction effect). The same LMM and GEE models were fitted, and the power was calculated as the percentage of simulations in which at least one null hypothesis was correctly rejected at the Bonferroni-corrected type I error level.

Simulation results

Table 3 presents the empirical type I error rates for LMM, standard GEE, and GEE with the Mancl & DeRouen small-sample correction (GEE-MD). The evaluation was conducted across various sample sizes ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} ), intraclass correlation coefficients (ICC), and the three models we reviewed. The LMM has slightly inflated type I error rates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\alpha )$$\end{document} but close to nominal level for most conditions.. For Model 1 ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} = 0.025), error rates were stable and ranged from 0.023 to 0.031. For the more complex Models 2 and 3 ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} ≈ 0.0167), the rates were appropriately lower, ranging from 0.016 to 0.020. The performance of the LMM was robust to changes in both sample size and ICC. The standard GEE approach showed inflation of the type I error rate, particularly at the smallest sample size ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30). For example, under Model 1 with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30, the GEE error rate was as high as 0.047 (nominal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} = 0.025), nearly double the target level. This liberal tendency was mitigated as the sample size increased, demonstrating the asymptotic nature of the method. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 120, the error rates, while still slightly elevated (e.g., 0.031), approached the nominal level. The application of the MD correction was highly effective in controlling the type I error of the GEE model. The corrected GEE produced error rates that were very similar to the LMM and were consistently close to the nominal alpha. For the condition where the standard GEE was most liberal (Model 1, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30), the MD corrected GEE had a much-improved error rate of 0.028. This demonstrates its utility in small-sample contexts. Table 3. Type I error (Bonferroni adjusted) of all models fitted with both linear mixed effect model (LMM) and generating estimating equations (GEE) under various intra-cluster correlation coefficientsICCDetailsLinear mixed model (LMM)Generalized estimating equations (GEE)GEE (Mancl & DeRouen)modelN \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ICC = 0.051300.0310.0260.0470.0420.0280.0231600.0280.0300.0390.0390.0280.03111200.0290.0230.0310.0280.0270.0232300.0200.0190.0180.0350.0340.0290.0170.0180.0152600.0170.0200.0140.0250.0280.0200.0170.0180.01421200.0160.0160.0170.0210.0210.0190.0170.0170.0163300.0200.0190.0190.0350.0340.0360.0170.0180.0193600.0170.0200.0190.0250.0280.0260.0170.0180.02031200.0160.0160.0200.0210.0210.0230.0170.0170.019ICC = 0.101300.0300.0270.0470.0420.0280.0231600.0280.0300.0360.0380.0260.02911200.0280.0250.0320.0280.0270.0252300.0200.0190.0180.0350.0340.0290.0180.0180.0152600.0180.0190.0140.0260.0250.0200.0180.0190.01421200.0170.0170.0170.0210.0200.0190.0160.0170.0163300.0200.0190.0220.0350.0340.0350.0180.0180.0193600.0180.0190.0190.0260.0250.0260.0180.0190.02031200.0170.0170.0180.0210.0200.0190.0160.0170.017ICC = 0.301300.0280.0240.0430.0380.0250.0211600.0260.0290.0330.0370.0240.02911200.0270.0250.0300.0290.0260.0262300.0200.0200.0180.0340.0350.0290.0160.0180.0152600.0200.0180.0140.0260.0240.0200.0180.0170.01421200.0170.0160.0170.0200.0200.0190.0160.0170.0163300.0200.0200.0210.0340.0350.0320.0160.0180.0173600.0200.0180.0180.0260.0240.0260.0180.0170.01831200.0170.0160.0170.0200.0200.0170.0160.0170.014N total number of participants, ICC intra-cluster correlation coefficient

When there is no interaction effect between the interventions, all models lead to unbiased estimates. In Table 4, as expected, across all analytical methods, power increased substantially with larger sample sizes and decreased as the strength of the intraclass correlation increased. A crucial trade-off between type I error control and statistical power was observed when comparing the three methods. The standard GEE model consistently produced the highest numerical power, especially at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30. For instance, in Model 2 with ICC = 0.30 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30, the GEE had a power of 0.511, whereas the LMM and MD-corrected GEE had powers of 0.454 and 0.393, respectively. However, this is a direct consequence of the inflated type I error rate. The LMM and GEE with MD correction, which both demonstrated control of type I error, yielded more conservative (and more accurate) power estimates. The GEE with the MD correction often had the lowest power. As another example, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} =30 and ICC = 0.05 (Model 1), the power was 0.740, compared to the LMM’s 0.792. This reflects the statistical “cost” of bias-corrected robust variance estimation in small-sample settings (Table 4). Table 4. Simulated power of detecting the treatment effects among all three models when there is no interaction effectICCDetailsLinear mixed model (LMM)Generalized estimating equations (GEE)GEE (Mancl & DeRouen)modelN \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ICC = 0.051300.7920.7900.8180.8180.7400.7411600.9780.9810.9790.9820.9730.97611201.0001.0001.0001.0001.0001.0002300.5130.5080.0180.5680.5490.0290.4400.4320.0152600.8400.8440.0140.8450.8510.0200.8130.8130.01421200.9940.9930.0170.9940.9930.0190.9930.9910.0163300.5130.5080.9720.5680.5490.9760.4400.4320.9523600.8400.8441.0000.8450.8511.0000.8130.8131.00031200.9940.9931.0000.9940.9931.0000.9930.9911.000ICC = 0.101300.7750.7730.7950.7960.7170.7191600.9710.9740.9740.9770.9660.97011201.0001.0001.0001.0001.0001.0002300.4940.4910.0180.5520.5370.0290.4300.4190.0152600.8270.8300.0140.8370.8400.0200.7950.7960.01421200.9920.9910.0170.9930.9910.0190.9920.9890.0163300.4940.4910.9620.5520.5370.9650.4300.4190.9373600.8270.8301.0000.8370.8401.0000.7950.7961.00031200.9920.9911.0000.9930.9911.0000.9920.9891.000ICC = 0.301300.7130.7140.7410.7370.6570.6591600.9520.9540.9550.9570.9430.94511200.9990.9990.9990.9990.9990.9992300.4540.4520.0180.5110.5030.0290.3930.3780.0152600.7840.7860.0140.7980.8020.0200.7540.7540.01421200.9860.9840.0170.9870.9850.0190.9840.9800.0163300.4540.4520.9280.5110.5030.9360.3930.3780.8953600.7840.7860.9980.7980.8020.9980.7540.7540.99831200.9860.9841.0000.9870.9851.0000.9840.9801.000N total number of participants, ICC intra-cluster correlation coefficient

Table 5 presents the estimated treatment effects from all models in cases where there is an interaction effect, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document} = 0.40. In this setting, Model 1 incorrectly omits the interaction term and therefore produces substantially biased estimates for the main effects of treatments A and B. Models 2 and 3, which correctly account for the interaction, yield unbiased estimates for all parameters. These findings were consistent regardless of sample size, estimation methods (GEE vs. LMM) or ICC, underscoring that the primary driver of accurate estimation is the correct specification of the model structure. Table 5. Estimates for parameters of interest when there is interaction between intervention A and BICCDetailsLinear mixed model (LMM)Generalized estimating equations (GEE)modelN \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ICC = 0.051301.0051.0011.0051.0011601.0001.0001.0001.00011200.9971.0010.9971.0012300.8070.8030.3950.8070.8030.3952600.7990.7990.4020.7990.7990.40221200.7990.8020.3960.7990.8020.3963300.8070.8032.0050.8070.8032.0053600.7990.7992.0010.7990.7992.00131200.7990.8021.9980.7990.8021.998ICC = 0.101301.0051.0011.0051.0011601.0001.0011.0001.00111200.9971.0010.9971.0012300.8080.8040.3950.8080.8040.3952600.7990.8000.4020.7990.8000.40221200.7990.8020.3960.7990.8020.3963300.8080.8042.0060.8080.8042.0063600.7990.8002.0010.7990.8002.00131200.7990.8021.9980.7990.8021.998ICC = 0.301301.0071.0021.0071.0021601.0011.0011.0011.00111200.9981.0010.9981.0012300.8090.8050.3950.8090.8050.3952600.7990.8000.4020.7990.8000.40221200.8000.8030.3960.8000.8030.3963300.8090.8052.0090.8090.8052.0083600.7990.8002.0020.7990.8002.00231200.8000.8031.9980.8000.8031.998N total number of participants, ICC intra-cluster correlation coefficient

The statistical power to detect the interaction effect was low across all conditions, especially in small samples (Table 6). With a sample size of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30, power was below 0.11 for all methods. Power increased with sample size but remained modest at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} =120, reaching about 0.32 to 0.33 for the LMM and GEE methods. This result highlights a well-known challenge. Studies often lack sufficient power to detect true, complex, interaction effects. A critical finding emerges when comparing the power for main effects between the Model 1 (which omits the interaction term) and the correctly specified Model 2. The Model 1 produced increased power for detecting main effects (e.g., 0.937 vs. 0.513 when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} = 30). This inflation occurs because when a true interaction is ignored, its variance is improperly absorbed by the main effect terms in the model. This makes the main effects artificially large and more statistically significant, leading to an incorrectly high power and potentially erroneous scientific conclusions. It creates an illusion of a strong, simple main effect when the reality is more nuanced. Table 6. Simulated power of detecting the treatment effects among all three models when there is an interaction effectICCDetailsLinear mixed model (LMM)Generalized estimating equations (GEE)GEE (Mancl & DeRouen)modelN \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{A}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{B}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{I}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{C}$$\end{document} ICC = 0.051300.9370.9410.9410.9430.9080.9071600.9990.9990.9991.0000.9990.99911201.0001.0001.0001.0001.0001.0002300.5130.5080.0760.5680.5490.1070.4400.4320.0572600.8400.8440.1510.8450.8510.1700.8130.8130.13321200.9940.9930.3170.9940.9930.3310.9930.9910.3033300.5130.5080.9990.5680.5490.9980.4400.4320.9943600.8400.8441.0000.8450.8511.0000.8130.8131.00031200.9940.9931.0000.9940.9931.0000.9930.9911.000ICC = 0.101300.9270.9280.9310.9320.8930.8921600.9990.9990.9990.9990.9980.99811201.0001.0001.0001.0001.0001.0002300.4940.4910.0750.5520.5370.1070.4300.4190.0572600.8270.8300.1510.8370.8400.1690.7950.7960.13321200.9920.9910.3170.9930.9910.3320.9920.9890.3033300.4940.4910.9980.5520.5370.9960.4300.4190.9933600.8270.8301.0000.8370.8401.0000.7950.7961.00031200.9920.9911.0000.9930.9911.0000.9920.9891.000ICC = 0.301300.8900.8920.9000.9000.8510.8551600.9950.9970.9950.9960.9940.99511201.0001.0001.0001.0001.0001.0002300.4540.4520.0760.5110.5030.1080.3930.3780.0562600.7840.7860.1500.7980.8020.1690.7540.7540.13321200.9860.9840.3170.9870.9850.3320.9840.9800.3033300.4540.4520.9890.5110.5030.9900.3930.3780.9813600.7840.7861.0000.7980.8021.0000.7540.7541.00031200.9860.9841.0000.9870.9851.0000.9840.9801.000N total number of participants, ICC intra-cluster correlation coefficient

Discussion

This article proposed and evaluated a multiple-baseline factorial design that can allow participants to receive all trial interventions, improving retention and reducing ethical concerns, while minimizing sample size requirements. Through simulations, we compared LMM, GEE, and GEE-MD for analyzing data from this design. We demonstrated that Model 1, which omits the interaction term, provides biased main effect estimates when a true interaction is present. As such, while Model 1 provides efficiency gains, we do not recommend its use due to the high risk of misspecification. Either Model 2 or 3 should be used to avoid biased conclusions, unless interactions are implausible. Critically, when performing sample size calculations, a model that accounts for potential interactions should be used to ensure the study is adequately powered for the more realistic scenario.

The primary limitation of the MBFD (for Model 1 and 2) is the potential for carryover or sequence effects, since participants receive a combined treatment (A and B) only after receiving an individual one (A or B). Model 3 attempts to account for this by explicitly modeling main effects, interactions, and the sequential structure. While this approach reduces bias relative to simpler models, it cannot fully eliminate the confounding introduced by the design (e.g., one can only have the averaged effects of AB and BA if sequence effects exist). Therefore, the parameters need to be interpreted based on the assumptions. A standard parallel-group factorial design would avoid sequencing concerns and provide the cleanest, unbiased estimates of both main and interaction effects. Therefore, the proposed MBFD is a pragmatic alternative in specific contexts, such as when participant recruitment is severely limited or when it is important for all individuals to receive all interventions over the course of the study. Whenever it is possible, we would still recommend a standard factorial design.

Our simulations comparing LMM and GEE approaches align with previous research (Westgate & Burchett, 2016). Standard GEE tend to have inflated type I error rates in small sample sizes. In contrast, GEE with the Mancl and DeRouen small-sample correction successfully maintained the nominal error rate. This highlights the necessity of using a small-sample correction for GEE in MFBD. While LMM provided slightly higher statistical power than the GEE-MD approach, the random effects structure must be correctly specified. Given that GEE is robust to the misspecification of the correlation structures, we recommend using GEE with a proven small-sample correction to ensure valid statistical inferences, especially when the number of individuals is limited.

This study has several limitations. First, our simulation parameters, including the ICC and the magnitude of treatment effects, were chosen for illustrative purposes and were not anchored to specific empirical data. Future applications of this design should use parameters and effect sizes justified by the literature. Second, our investigation of statistical methods was not exhaustive. We focused on the MD correction for GEE, but other methods like the Kenward-Roger correction for LMM or the Kauermann and Carroll correction for GEE warrant investigation (Kauermann & Carroll, 2001; Kenward & Roger, 1997). More recent research showed that a robust variance estimator can work well under LMM for longitudinal and clustered data (Ouyang et al., 2024). Finally, future work could expand the framework to directly compare the combined intervention to individual interventions (e.g., AB vs. A), which is often a key question of interest in factorial trials.

Conclusion

In this study, we described a novel individually randomized MBFD to investigate two or more interventions against a control or SoC when the sample size is limited. Given the sequential exposure to treatments and the potential confounding from treatment order and carryover, we carefully considered the causal assumptions required to interpret effects, particularly when estimating combined treatment effects. We recommend using this design when the expected recruitment sample size is small and a conventional factorial design cannot provide sufficient power. The estimand of interest must be clarified, and either LMM or GEE, with a small sample correlation, is the most suitable analysis method.