# Learn-As-you-GO (LAGO) Trials: Optimizing Treatments and Preventing Trial
Failure Through Ongoing Learning

Ante<EMAIL_ADDRESS>[ Donna
<EMAIL_ADDRESS>[ Daniel
<EMAIL_ADDRESS>0000-0002-9770-827X [ Judith J.
<EMAIL_ADDRESS>[ Department of Mathematics and
Statistics, Boston Universitypresep=, ]e1,e4  Department of Biostatistics,
Yale Universitypresep=, ]e2  Department of Statistics and Operations Research,
Tel Aviv Universitypresep=, ]e3

###### Abstract

It is well known that changing the intervention package while a trial is
ongoing does not lead to valid inference using standard statistical methods.
However, it is often necessary to adapt, tailor, or tweak a complex
intervention package in public health implementation trials, especially when
the intervention package does not have the desired effect. This article
presents conditions under which the resulting analyses remain valid even when
the intervention package is adapted while a trial is ongoing. Our results on
such Learn-As-you-GO (LAGO) studies extend the theory of LAGO for binary
outcomes following a logistic regression model (Nevo, Lok and Spiegelman,
2021) to LAGO for continuous outcomes under flexible conditional mean model.
We derive point and interval estimators of the intervention effects and ensure
the validity of hypothesis tests for an overall intervention effect. We
develop a confidence set for the optimal intervention package, which achieves
a pre-specified mean outcome while minimizing cost, and confidence bands for
the mean outcome under all intervention package compositions. This work will
be useful for the design and analysis of large-scale intervention trials where
the intervention package is adapted, tailored, or tweaked while the trial is
ongoing.

62K99,

62L99,

62F05,

62F10,

62F12,

62J12,

Adaptive trial designs,

dependent sample,

large-scale intervention trials,

public health,

implementation trials,

###### keywords:

[class=MSC]

###### keywords:

, , and

##  1 Introduction

Traditionally regulators have not allowed adaptive trial designs to be used
for evidence in the drug approval process, but currently, the Food and Drug
Administration (FDA) has put forward guidelines for adaptive clinical trials
(FDA, 2016, 2019) and encourages trialists interested in adaptive trials to
contact the FDA in the planning stage (FDA, 2019). While current adaptive
clinical trial designs allow for changing randomization probabilities and
dropping treatment arms (FDA, 2019), they do not allow for changes to an
intervention package composition based on outcomes collected during the
earlier stages of a trial. This limitation may have contributed to cases where
large-scale intervention trials have “failed” (Semrau et al., 2017; Fogel,
2018; Stensland et al., 2014).

Learn-As-you-GO (LAGO) studies consist of $K>1$ stages. As data are being
collected, the later-stage intervention package is systematically “learned”
based on results from previous stages. Nevo, Lok and Spiegelman (2021)
provided a methodology for LAGO studies with binary outcomes. They assumed
that a logistic regression model holds for the probability of success given
the intervention package components, for which they prove consistency and
asymptotic normality. In LAGO trials, learning assumes that as the number of
patients in prior stages gets large, the recommended intervention in the next
stage converges in probability to a fixed intervention. They also demonstrated
the validity of hypothesis tests for the overall intervention effect.
Moreover, they proved that the optimal intervention, $\boldsymbol{x}^{opt}$,
which attains a pre-specified outcome probability while minimizing cost, can
be “learned”, in the sense that we can find $\hat{\boldsymbol{x}}^{opt}$ with
$\hat{\boldsymbol{x}}^{opt}\stackrel{{\scriptstyle
P}}{{\rightarrow}}\boldsymbol{x}^{opt}$ as the number of patients in each
stage gets large.

An adaptive clinical design popular in implementation science (Collins et al.,
2011) that can be compared to LAGO is the multiphase optimization strategy
(MOST) (Collins, Murphy and Strecher, 2007; Collins, Nahum-Shani and Almirall,
2014). Although both LAGO and MOST aim to identify the optimal intervention
package and evaluate its impact, the designs are different. MOST consists of
three phases: preparation, optimization, and evaluation. In the optimization
phase, the optimal intervention package is determined in a short-term
factorial design usually of intermediate outcomes, and in the evaluation
phase, the effect of the previously identified intervention package is
independently assessed through a randomized controlled trial (RCT), the
analysis of which does not use the outcomes from the optimization phase or any
occurrences subsequently. In contrast, in LAGO trials, the composition of the
intervention package is updated at the end of each stage and is dependent on
the outcomes from previous stages. LAGO enables researchers to systematically
modify the intervention package while the trial is ongoing while preserving
Type I error rate and pre-specified power. Since all outcomes are used in the
final analysis, as further discussed in Section 7, LAGO can achieve higher
power compared to MOST when the same strategy is used to identify the optimal
intervention package.

LAGO has been identified as a promising approach in the field of
implementation science (Beidas et al., 2022). According to these authors, LAGO
can improve the alignment between implementation strategies and partner needs
and contexts. Additionally, the user-centered design and approach of LAGO
enables the optimization of implementation strategies and calibration of
implementation support based on demonstrated need.

This article extends LAGO to continuous outcomes. To achieve this, we adopt a
Generalized Linear Model (GLM) framework. The GLM framework in this context
refers to any semi-parametric model for the conditional mean of a random
outcome on other covaraites that is linear in a link function, without any
additional restrictions on the outcome distribution. When analyzing a LAGO
study, we cannot condition on the later-stage interventions since that would
imply conditioning on a function of the outcomes from earlier stages. In
addition, variations in the intervention components is needed for LAGO to
identify the treatment effect parameters. Extending LAGO to continuous
outcomes is not straightforward, as Nevo, Lok and Spiegelman (2021) use a
novel coupling argument to prove the asymptotic properties of their
estimators, and such coupling approach cannot be generalized to continuous
outcomes. To overcome this challenge, we assume throughout this article that
the errors in the GLM (McCullagh and Nelder, 2019) are independent of the
composition of the intervention package. Section 4 outlines the comprehensive
steps implemented to surmount the difficulty in generalizing LAGO to
continuous outcomes. Details are provided in Appendix B.

The aims of a LAGO study with continuous outcomes are 1$)$ to estimate the
impact of the individual intervention package components on the outcome mean,
2$)$ to test for an overall intervention effect when the intervention package
has been “learned,” and 3$)$ to estimate which intervention package,
$\boldsymbol{x}^{opt}$, will minimize the cost while yielding a pre-specified
mean outcome. We solve these aims for various GLMs with independent,
identically distributed errors.

Section 8 applies LAGO to the BetterBirth study (Hirschhorn et al., 2015;
Semrau et al., 2017). The BetterBirth study was a costly failed trial, led by
Harvard researcher Atul Gawande, which aimed to improve maternal and child
health around the time of birth. By mimicking applying LAGO, we identify the
optimal intervention package that can increase the percentage of essential
birth practices (EBPs) performed to a pre-specified target goal while
minimizing cost. Our analysis differs from that in Nevo, Lok and Spiegelman
(2021) because they focused on a single binary outcome, oxytocin administered
immediately after delivery, as opposed to the continuous outcome of the
percentage of EBPs performed. Furthermore, since we consider the percentage of
EBPs as a continuous outcome, this represents a unique application of our LAGO
method that could not have been accomplished using the existing LAGO method.

This article is organized as follows. Section 2 details the setting and
notation. Section 3 describes the estimating equations for the proposed
estimator for a LAGO study with a GLM. Section 4 outlines the proofs of the
asymptotic properties of the proposed estimator. Details are presented in
Appendix B. Section 5 describes how to test for an overall intervention effect
in a LAGO study. Section 6 describes the confidence sets and confidence bands
for the optimal intervention. Section 7 describes simulations. Section 8
describes the LAGO analysis of the BetterBirth study. Section 9 discusses our
findings and future research topics.

##  2 Setting, notation, and assumptions

Let $Y_{ij}$ be the continuous outcome for patient $i$ in center $j$. The
optimal intervention is defined as the intervention that results in the mean
outcome reaching a pre-specified goal $\theta$, while minimizing cost. Suppose
the multi-component intervention package $\boldsymbol{x}$ has $p$ components.
Let $C(\boldsymbol{x})$ be a known cost function, in dollars, of the
intervention package $\boldsymbol{x}$. In healthcare settings, a linear cost
function is often reasonable to assume. $C(\boldsymbol{x})$ will then include
a fixed cost plus the sum of the product between dosage and the unit cost for
each of the intervention components. We also consider a cubic cost function,
which allows for an initial economy of scale followed by increased costs when
the component levels exceed a threshold (see Appendix E.5 for details). This
type of cubic cost function is particularly applicable in healthcare settings,
where the marginal cost of a product may rise prohibitively as local supplies
are depleted. In addition, in the BetterBirth study, one of the intervention
components is the duration of the on-site intervention launch, measured in
days (see Section 8 for more details on the BetterBirth study). Here, a cubic
cost function is suitable because extending the launch beyond three days would
disrupt care at the participating healthcare centers, making it infeasible.
Typically, less feasible doses are associated with higher costs.

Let $\boldsymbol{z}_{j}$ be fixed center-specific characteristics, such as
hospital district or hospital birth volume, that may be related to the outcome
of interest.

######  Assumption 1.

We assume that the expected outcome of an individual $i$ in center $j$ with
characteristic $\boldsymbol{z}_{j}$ under recommended intervention package
$\boldsymbol{X}_{j}=\boldsymbol{x}_{j}$ and actual intervention package
$\boldsymbol{A}_{j}=\boldsymbol{a}_{j}$,
$E\left(Y_{ij}|\boldsymbol{A}_{j}=\boldsymbol{a}_{j},\boldsymbol{X}_{j}=\boldsymbol{x}_{j},\boldsymbol{z}_{j};\boldsymbol{\beta}\right)$
only depends on the recommended intervention $\boldsymbol{x}_{j}$ through the
actual intervention $\boldsymbol{a}_{j}$ and follows a GLM

$g\left(E\left(Y_{ij}|\boldsymbol{A}_{j}=\boldsymbol{a}_{j},\boldsymbol{X}_{j}=\boldsymbol{x}_{j},\boldsymbol{z}_{j};\boldsymbol{\beta}\right)\right)=\beta_{0}+\boldsymbol{\beta}_{1}^{T}\boldsymbol{a}_{j}+\boldsymbol{\beta}_{2}^{T}\boldsymbol{z}_{j},$
(1)

where $g()$ is a twice continuously differentiable link function.
$\boldsymbol{\beta}^{T}=\left({\beta}_{0},\boldsymbol{\beta}_{1}^{T},\boldsymbol{\beta}_{2}^{T}\right)$
are unknown parameters to be estimated from the data. Typically, our main
interest is in the vector $\boldsymbol{\beta}_{1}$, and the null hypothesis of
no intervention effect, $H_{0}:\boldsymbol{\beta}_{1}=\boldsymbol{0}$.

The optimal intervention package for center $j$ with baseline covariates
$\boldsymbol{z}_{j}$ is the solution to the center-specific optimization
problem

$\operatorname{Min}_{\boldsymbol{x}_{j}}C\left(\boldsymbol{x}_{j}\right)\;\;\text{
subject to
}\;\;E\left(Y_{ij}|\boldsymbol{x}_{j},\boldsymbol{z}_{j}\;;\boldsymbol{\beta}\right)\geq\theta,$
(2)

with each component $p$ of $\boldsymbol{x}_{j}$ in a pre-determined interval
$\left[L_{p},U_{p}\right]$, where $p=1,...,P$. Alternatively, instead of an
absolute goal, the optimal intervention package could e.g. aim to increase the
mean outcome for each center $j$ by a pre-specified goal $\Delta\theta$. In
what follows, we focus on case (2). In many settings, the intervention package
implemented will not be center-specific, in which case all $\boldsymbol{z}$
are ignored (or set to 0). We assume that there is a unique solution to
equation (2). If we consider a linear cost function, and the intervention
consists of two components with marginal costs $c_{1}$ and $c_{2}$, then if
$\beta_{11}/c_{1}\neq\beta_{12}/c_{2}$, the solution to equation (2) is
unique. Part of the aims of a LAGO study is to identify the solution,
$\boldsymbol{x}^{opt}_{j}$, that solves equation (2).

For simplicity, we present the theory for LAGO studies with $K=2$ stages.
Appendix D extends the theory to $K>2$ stages.

In each stage $k$, $k=1$ or $2$, $n_{j}^{(k)}$ patients in center $j$,
$j=1,\ldots,J^{(k)}$ are enrolled, with $J^{(k)}$ fixed. In a randomized
controlled LAGO trial (c-LAGO), centers may be randomized to either the
intervention or control arm. In a quasi-experimental LAGO study, pre-study
outcome data might be collected (before-after or ba-LAGO). Let
$n^{(k)}=\sum_{j=1}^{J^{(k)}}n_{j}^{(k)}$ be the number of patients in stage
$k$, $n=n^{(1)}+n^{(2)}$ the total number of patients across the 2 stages, and
$\alpha_{jk}=\lim_{n\rightarrow\infty}n_{j}^{(k)}/n\;>0$ , which assumes that
the ratio between the sample size for center $j$ in stage $k$ and the total
sample size converges to a non-zero constant as $n$ goes to infinity. This
assumption is a reasonable approximation if the number of patients in each
center $j$ of stages 1 and 2 is large, as in most large scale implementation
trials. At stage $1$, the initial intervention package for each center,
$\boldsymbol{x}_{j}^{(1)}$, is recommended by the study investigators’ best
guess, and/or based on pilot data before the study starts. In stage $2$, the
recommended intervention package $\boldsymbol{X}_{j}^{(2,n^{(1)})}$ for center
$j$ is estimated from stage $1$ data and other accumulated knowledge, possibly
including summary measures of e.g. qualitative information from providers,
patients and other stakeholders. The superscript $(2,n^{(1)})$ indicates that
$\boldsymbol{X}_{j}^{(2,n^{(1)})}$ depends on the data from the $n^{(1)}$
patients in stage 1. Often, $\boldsymbol{X}_{j}^{(2,n^{(1)})}$ will solve
equation (2) with the estimator $\hat{\boldsymbol{\beta}}^{(1)}$ for
$\boldsymbol{\beta}$ based on the stage 1 data (see Remark 1). Let
$\boldsymbol{A}_{j}^{(2,n^{(1)})}$ be the actual intervention package
implemented in center $j$ in stage $2$. Under perfect adherence,
$\boldsymbol{A}_{j}^{(2,n^{(1)})}=\boldsymbol{X}_{j}^{(2,n^{(1)})}$. In
practice, both stage 1 and stage 2 centers may not adhere completely to all
the intervention components (unplanned variation or uv-LAGO). We assume that
to the extent it is not random, adherence depends on the recommendation
$\boldsymbol{X}_{j}^{(2,n^{(1)})}$ and on the stage $2$ center characteristics
$\boldsymbol{z}_{j}^{(2)}$, but not further on other predictors of the
outcome.

######  Remark 1.

One approach to determining recommended interventions for stage 2.
One way to determine the recommended interventions
$\boldsymbol{X}_{j}^{(2,n^{(1)})}$ for stage 2, denoted by
$\hat{\boldsymbol{x}}^{opt,(2,n^{(1)})}_{j}$, is by using a function $f$. This
function takes stage 1 outcomes, the pre-specified goal $\theta$, and stage 2
center-specific characteristics $\boldsymbol{z}_{j}^{(2)}$ as input, solving
the optimization problem of equation (2) using the stage 1-based estimate
$\hat{\boldsymbol{\beta}}^{(1)}$ in place of $\boldsymbol{\beta}$, and returns
center-specific recommended interventions as output. Therefore,
$\hat{\boldsymbol{x}}^{opt,(2,n^{(1)})}_{j}=f(\hat{\boldsymbol{\beta}}^{(1)};\boldsymbol{z}_{j}^{(2)})$
and $f$ solves equation (2).

The optimization algorithm used to solve Equation (2) varies depended upon
form of the cost function, $C(\boldsymbol{x})$. For linear cost functions, we
first rank the intervention components based on their cost efficiency,
$\hat{\beta}_{1p}/c_{p}$, where $c_{p}$ is the unit cost for the $pth$
intervention component. Suppose that the intervention consists of two
components, $p\in\\{1,2\\}$. We increase the more cost-efficient component,
while keeping the other component at its minimum, to check whether the
objective $\theta$ as described in Equation (2) can be achieved. If the
objective cannot be met under this setting, we set the more cost-efficient
component to its maximum value and increase the less cost-efficient component
until the objective $\theta$ as described in Equation (2) is achieved. If the
target value $\tilde{\theta}$ cannot be achieved, a function
$g(\hat{\boldsymbol{\beta}})$ is employed to determine the recommended
intervention. This function ensures that the recommended intervention is
continuous with respect to each $\hat{\beta}_{1p}$. Nevo, Lok and Spiegelman
(2021) described this first, and further details can be found in Section 5.1
of their Supplementary Material.

In the case of cubic cost functions, we construct a grid with small increments
(e.g., 0.01) for each component. We then search for the combination of the two
components that satisfies the target criterion $\theta$ as outlined in
Equation (2), while minimizing the cubic cost function.

For clarity of the expositions, below we will work with
$\hat{\boldsymbol{x}}^{opt,(2,n^{(1)})}_{j}$. Let $Y_{ij}^{(1)}$ be the
outcome of participant $i$ in center $j$ of stage 1. Let
$\boldsymbol{a}_{j}^{(1)}=h_{j}^{(1)}\left(\boldsymbol{x}_{j}^{(1)}\right)$ be
the actual intervention package for center $j$ of stage 1, where $h_{j}^{(1)}$
is a continuous, deterministic function for each center $j$ in stage 1. Let
$\boldsymbol{Y}_{j}^{(1)}=\bigl{(}Y_{1j}^{(1)},\ldots,Y_{n_{j}^{(1)}j}^{(1)}\bigl{)}$
be the outcomes of patients $1,\ldots,n_{j}^{(1)}$ in center $j$ of stage 1.
Let
$\overline{\boldsymbol{a}}^{(1)}=\left(\boldsymbol{a}_{1}^{(1)},\ldots,\boldsymbol{a}_{J^{(1)}}^{(1)}\right)$,
$\overline{\boldsymbol{z}}^{(1)}=\left(\boldsymbol{z}_{1}^{(1)},\ldots,\boldsymbol{z}_{J^{(1)}}^{(1)}\right)$,
and
$\overline{\boldsymbol{Y}}^{(1)}=\left(\boldsymbol{Y}_{1}^{(1)},\ldots,\boldsymbol{Y}_{J^{(1)}}^{(1)}\right)$
be the actual interventions, center-specific characteristics and outcomes for
each center $1,\ldots,J^{(1)}$ of stage 1, respectively. Additionally, let
$\overline{\hat{\boldsymbol{x}}}^{opt,\left(2,n^{(1)}\right)}=\Bigl{(}\hat{\boldsymbol{x}}_{1}^{opt,\left(2,n^{(1)}\right)},\ldots,\hat{\boldsymbol{x}}_{J^{(2)}}^{opt,\left(2,n^{(1)}\right)}\Bigl{)}$
be the recommended interventions for center $1,\ldots,J^{(2)}$ of stage 2.

Let $Y_{ij}^{(2,n^{(1)})}$, $i=1,\cdots,n_{j}^{(2)}$ be the outcome of
participant $i$ in center $j$ of stage 2. Let
$\boldsymbol{Y}_{j}^{(2,n^{(1)})}=\Bigl{(}Y_{1j}^{(2,n^{(1)})},\ldots,Y_{n_{j}^{(2)}j}^{(2,n^{(1)})}\Bigl{)}$
be the outcomes of patients $1,\ldots,n_{j}^{(2)}$ in center $j$ of stage 2.
Let
$\boldsymbol{A}_{j}^{(2,n^{(1)})}=h_{j}^{(2)}\left(\hat{\boldsymbol{x}}^{opt,(2,n^{(1)})}_{j}\right)$
be the actual intervention package for center $j$ of stage 2, where
$h_{j}^{(2)}$ is a continuous, deterministic function for each center $j$ in
stage 2. Let
$\overline{\boldsymbol{A}}^{\left(2,n^{(1)}\right)}=\Bigl{(}\boldsymbol{A}_{1}^{\left(2,n^{(1)}\right)},\ldots,\boldsymbol{A}_{J^{(2)}}^{\left(2,n^{(1)}\right)}\Bigl{)}$,
$\overline{\boldsymbol{z}}^{(2)}=\left(\boldsymbol{z}_{1}^{(2)},\ldots,\boldsymbol{z}_{J^{(2)}}^{(2)}\right)$,
and
$\overline{\boldsymbol{Y}}^{\left(2,n^{(1)}\right)}=\Bigl{(}\boldsymbol{Y}_{1}^{\left(2,n^{(1)}\right)},\ldots,\boldsymbol{Y}_{J^{(2)}}^{\left(2,n^{(1)}\right)}\Bigl{)}$
be the actual interventions, center-specific characteristics and outcomes for
each center $1,\dots,J^{(2)}$ of stage 2, respectively.

######  Assumption 2.

Conditionally on
$\overline{\hat{\boldsymbol{x}}}^{opt,\left(2,n^{(1)}\right)}$,
$\bigl{(}\overline{\boldsymbol{A}}^{\left(2,n^{(1)}\right)},\overline{\boldsymbol{Y}}^{\left(2,n^{(1)}\right)}\bigl{)}$
is independent of the stage 1 data
$\bigl{(}\overline{\boldsymbol{a}}^{(1)},\overline{\boldsymbol{Y}}^{(1)}\bigl{)}$.
That is, learning from data from earlier stages is only through the
determination of the recommended intervention.

######  Assumption 3.

For each center $j=1,\ldots,J^{(2)}$, the stage 2 recommended intervention
$\hat{\boldsymbol{x}}^{opt,(2,n^{(1)})}_{j}$ converges in probability to a
center-specific limit $\boldsymbol{x}_{j}^{(2)}$.

######  Remark 2.

Assumption 3 holds e.g. when the recommended intervention is a continuous
function of either 1) the maximum likelihood estimator (MLE) for
$\boldsymbol{\beta}$; 2) a solution to generalized estimating equation (GEE,
Liang and Zeger (1986)); 3) averages based on all stage 1 patients. Thus,
Assumption 3 holds if $\hat{\boldsymbol{x}}^{opt,(2,n^{(1)})}_{j}$ solves
equation (2) and the solution to (2) is unique for the true
$\boldsymbol{\beta}$.

Under Assumption 3, the definition of $h_{j}$ that maps recommended
interventions to actual interventions in center $j$, and the Continuous
Mapping Theorem, we conclude that
$\boldsymbol{A}_{j}^{(2,n^{(1)})}=h_{j}^{(2)}\bigl{(}\hat{\boldsymbol{x}}^{opt,(2,n^{(1)})}_{j}\bigl{)}$
converges in probability to
$\boldsymbol{a}_{j}^{(2)}=h_{j}^{(2)}\bigl{(}\boldsymbol{x}^{(2)}_{j}\bigl{)}$.

######  Assumption 4.

The covariates $\boldsymbol{z}$, the outcomes $Y$, and the parameter space for
$\boldsymbol{\beta}$ all take values in a compact space.

######  Assumption 5.

Let $\epsilon_{ij}$ be the error term for individual $i$ in center $j$ of the
GLM from Assumption 1. That is,
$Y_{ij}=g^{-1}\left(\beta_{0}+\boldsymbol{\beta}_{1}^{T}\boldsymbol{a}_{j}+\boldsymbol{\beta}_{2}^{T}\boldsymbol{z}_{j}\right)+\epsilon_{ij}$,
with $E\left(\epsilon_{ij}|\boldsymbol{a}_{j},\boldsymbol{z}_{j}\right)=0$.
The distribution of the errors $\epsilon_{ij}$ is independent of the
composition of the intervention package $\boldsymbol{a}_{j}$, although it may
depend on $\boldsymbol{z}_{j}$. We denote
$\sigma^{2}(\boldsymbol{z}_{j})=VAR(\epsilon_{ij}|\boldsymbol{z}_{j})$. Often,
the $\epsilon_{ij}$ will be assumed to be iid, with distribution not dependent
on the $\boldsymbol{z}_{j}$.

##  3 Estimating $\boldsymbol{\beta}$ and the optimal intervention package

This section describes the estimating equations for the proposed estimator
$\hat{\boldsymbol{\beta}}$, and how to subsequently use
$\hat{\boldsymbol{\beta}}$ to estimate the optimal intervention package. The
adaption of the intervention package based on prior outcomes causes dependence
between stages, invalidating standard statistical theory. The usual (non-LAGO)
estimator $\hat{\boldsymbol{\beta}}$ for the true parameter
$\boldsymbol{\beta}^{*}$ is the solution to the GEE under an independence
working correlation structure (Liang and Zeger, 1986)

$\displaystyle 0=\boldsymbol{U}^{(g)}(\boldsymbol{\beta})$ (3)
$\displaystyle=\frac{1}{n}\left\\{\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right.$
$\displaystyle\qquad+\left.\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(Y_{ij}^{\left(2,n^{(1)}\right)}-g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right\\}.$

We estimate $\boldsymbol{\beta}$ in a LAGO study by the solution to equation
(3). The superscript $(g)$ in $\boldsymbol{U}^{(g)}(\boldsymbol{\beta})$
reminds us that $\boldsymbol{U}^{(g)}(\boldsymbol{\beta})$ are estimating
equations for LAGO GLM with a general link function. Asymptotic theory for
$\hat{\boldsymbol{\beta}}$ is complicated by the fact that the stage 2
interventions $\overline{\boldsymbol{A}}^{\left(2,n^{(1)}\right)}$ depend on
the stage 1 outcomes $\overline{\boldsymbol{Y}}^{(1)}$, so the two terms in
(3) are not independent.

Section 4 shows that despite these dependencies, under the assumptions of
Section 2, the estimator $\hat{\boldsymbol{\beta}}$ solving equation (3) is
both consistent and asymptotically normal.

To identify the final recommended intervention package for any given center,
we use the function $f$ from Remark 1, along with the stage 1 and 2 outcome
data, the pre-specified goal $\theta$, and the center-specific characteristics
to solve the optimization problem given by (2), but with the final estimator
$\hat{\boldsymbol{\beta}}$ in place of $\boldsymbol{\beta}$. The recommended
intervention package for a center with any specified center characteristics
$\boldsymbol{z}$ will be returned as output. That is,
$\hat{\boldsymbol{x}}^{opt}=f(\hat{\boldsymbol{\beta}};\boldsymbol{z})$and $f$
solves equation (2).

##  4 Asymptotic Properties of $\hat{\boldsymbol{\beta}}$

We describe the asymptotic properties of the final estimator
$\hat{\boldsymbol{\beta}}$. Similar to the previous sections, consider a two-
stage LAGO design. Section 4 outlines the proofs of consistency and asymptotic
normality; Appendix B provides the detailed proofs.

######  Theorem 1.

(Consistency). Under Assumptions 1 – 5,
$\hat{\boldsymbol{\beta}}\xrightarrow{P}\boldsymbol{\beta}^{*}$.

Nevo, Lok and Spiegelman (2021) provide a proof of consistency for
$\hat{\boldsymbol{\beta}}$ under logistic regression. Here, the general link
function $g()$ does not need to be the canonical link function, and the
outcomes $Y_{ij}$ are continuous.

To prove consistency of $\hat{\boldsymbol{\beta}}$, we show that in spite of
the fact that equation (3) does not consist of i.i.d. terms, Theorem $5.9$ of
Van der Vaart (2000) can be used. We show that the two conditions of this
theorem are satisfied. For the first condition of Theorem $5.9$ of Van der
Vaart (2000), we show that

$\small{\sup_{\boldsymbol{\beta}}\|\boldsymbol{U}^{(g)}(\boldsymbol{\beta})-\boldsymbol{u}^{(g)}(\boldsymbol{\beta})\|\xrightarrow{P}0,}$
(4)

where $\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$ are the expected values of
the estimating equations under the limiting design, in particular,

$\displaystyle\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$
$\displaystyle=\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right]$
(5)
$\displaystyle+\sum_{j=1}^{J^{(2)}}\alpha_{j2}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})-g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right],$

where (see Section 2), $\alpha_{jk}=\lim_{n\rightarrow\infty}n_{j}^{(k)}/n$
and $\boldsymbol{A}_{j}^{(2,n^{(1)})}\xrightarrow{P}\boldsymbol{a}_{j}^{(2)}$.
(4).

To prove equation (4), Appendix B.1 shows that
$\boldsymbol{U}^{(g)}(\boldsymbol{\beta})-\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$
can be decomposed into five distinct components. Let $\epsilon_{ij}^{(2)}$ be
the error that patient $i$ in center $j$ would have experienced under the
limiting intervention $\boldsymbol{a}_{j}^{(2)}$. By Assumption 5, replacing
the $\epsilon_{ij}^{(2,n^{(1)})}$ by the error terms $\epsilon_{ij}^{(2)}$
does not change the distribution of the part in
$\boldsymbol{U}^{(g)}(\boldsymbol{\beta})-\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$
that includes $\epsilon_{ij}^{(2,n^{(1)})}$. Combining the error replacing
approach with the concept of Donsker classes from empirical process theory,
the supremum of $\boldsymbol{\beta}$ for each of the five components converges
to 0 in probability. Then by triangle inequality, equation (4) holds.

In order for LAGO to work properly, we need variations in the intervention
components to identify the treatment effect parameters. The uniqueness of
$\boldsymbol{\beta}^{*}$ as a maximizer or zero of
$\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$ (equation (5)) has been studied by
various authors, see e.g. Chapter 2.2 of Fahrmeir and Tutz (2013).
$\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$ is the same as if the interventions
were fixed before the study, so if there is enough variation in the
intervention components, the second condition in Theorem 5.9 of Van der Vaart
(2000) is also satisfied and we conclude that $\hat{\boldsymbol{\beta}}$ is
consistent.

######  Theorem 2.

(Asymptotic Normality). Under Assumptions 1 – 5,

$\small{\sqrt{n}\left(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}\right)\xrightarrow{D}N\left(0,J\left(\boldsymbol{\beta^{*}}\right)^{-1}V\left(\boldsymbol{\beta^{*}}\right)J\left(\boldsymbol{\beta^{*}}\right)^{-1}\right),}$
(6)

where

$\displaystyle J(\boldsymbol{\beta}^{*})$
$\displaystyle=\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}+\sum_{j=1}^{J^{(2)}}\alpha_{j2}\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2},$ $\displaystyle V\left(\boldsymbol{\beta}^{*}\right)=$
$\displaystyle\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\sigma^{2}(\boldsymbol{z}_{j}^{(1)})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j=1}^{J^{(2)}}\alpha_{j2}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2}\sigma^{2}(\boldsymbol{z}_{j}^{(2)}),$

and ⊗2 is the Kronecker product.

The variance in equation (6) can be estimated by replacing
$\boldsymbol{\beta}^{*}$, $\boldsymbol{a}_{j}^{(2)}$, $\alpha_{j1}$, and
$\alpha_{j2}$ with $\hat{\boldsymbol{\beta}}$,
$\boldsymbol{A}_{j}^{(2,n^{(1)})}$, ${n_{j}^{(1)}}/{n}$, and
${n_{j}^{(2)}}/{n}$, respectively. For estimating
$\sigma^{2}(\boldsymbol{z}_{j})$ and exact formulas of
$\hat{J}(\hat{\boldsymbol{\beta}})$ and $\hat{V}(\hat{\boldsymbol{\beta}})$,
see equations (21) and (22) in Appendix B.2.

Nevo, Lok and Spiegelman (2021) provide a proof of asymptotic normality for
$\hat{\boldsymbol{\beta}}$ under logistic regression; however, the novel
coupling arguments they used for binary outcomes cannot be generalized to
continuous outcomes. In LAGO trials, standard theory is not applicable as the
stage 2 intervention package depends on stage 1 outcomes. Specifically,
conditioning on the stage 2 intervention package means conditioning on
functions of stage 1 outcomes.

To address this issue, we assume that the errors in the GLM are independent of
the intervention package composition (see Assumption 5). The proof of Theorem
2 employs a combination of strategies and concepts: building on the
established Theorem 1, utilizing the Mean Value Theorem, employing the error-
replacement strategy from the proof of Theorem 1, invoking Donsker classes
from empirical process theory, and leveraging Lévy’s Continuity Theorem for
characteristic functions. Ultimately, Appendix B.2 demonstrates that
$\sqrt{n}\bigl{(}\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}\bigl{)}$ has
the same asymptotic distribution as the estimator for a fixed two-stage design
with interventions $\boldsymbol{a}_{j}^{(2)}$ decided on before the trial.

By the Mean Value Theorem (similar to Theorem 5.21 of Van der Vaart (2000)),

$\small{\sqrt{n}\left(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}\right)=-\sqrt{n}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}\left({\boldsymbol{\beta}}\right)\right)^{-1}\boldsymbol{U}^{(g)}\left(\boldsymbol{\beta}^{*}\right),}$

where for each row of
$\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}({\boldsymbol{\beta}})$,
$\tilde{\boldsymbol{\beta}}$ takes a possibly row-dependent value between
$\hat{\boldsymbol{\beta}}$ and $\boldsymbol{\beta}^{*}$. By Theorem 1, the
Continuous Mapping Theorem and further arguments, Appendix B.2 shows that
$\bigl{(}-\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}\left({\boldsymbol{\beta}}\right)\bigl{)}$
converges in probability to $J\left(\boldsymbol{\beta^{*}}\right)$. Next, we
derive

$\small{\sqrt{n}\
\boldsymbol{U}^{(g)}(\boldsymbol{\beta}^{*})=\boldsymbol{U}_{1,n}^{(g)}+\boldsymbol{U}_{2\\_1,n}^{(g)}+\boldsymbol{U}_{2\\_2,n}^{(g)},}$
(7)

where

$\displaystyle\boldsymbol{U}_{1,n}^{(g)}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\epsilon_{ij}^{(1)},$
$\displaystyle\boldsymbol{U}_{2\\_2,n}^{(g)}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\epsilon_{ij}^{(2,n^{(1)})},$

and the explicit form of $\boldsymbol{U}_{2\\_1,n}^{(g)}$ is given by equation
(29) in Appendix B.2.

Appendix B.2 shows that $\boldsymbol{U}_{2\\_1,n}^{(g)}\xrightarrow{P}0$. For
$\boldsymbol{U}_{2\\_2,n}^{(g)}$, similar to the proof of Theorem 1 in
Appendix B.1, by Assumption 5, replacing $\epsilon_{ij}^{(2,n^{(1)})}$ by the
error terms $\epsilon_{ij}^{(2)}$ under the limiting interventions
$\boldsymbol{a}_{j}^{(2)}$ does not change the distribution of
$\boldsymbol{U}_{2\\_2,n}^{(g)}$. We denote
${\boldsymbol{U}_{2\\_2,n}^{(g)\;*}}$ as the new term with error terms
$\epsilon_{ij}^{(2)}$ and interventions $\boldsymbol{a}_{j}^{(2)}$. This
implies that $\sqrt{n}\ \boldsymbol{U}^{(g)}(\boldsymbol{\beta}^{*})$ has the
same asymptotic distribution as
$\boldsymbol{U}_{1,n}^{(g)}+{\boldsymbol{U}_{2\\_2,n}^{(g)\;*}}$, which are
the estimating equations for a fixed two-stage design with
$\boldsymbol{a}_{1}^{(1)},\cdots,\boldsymbol{a}_{J^{(1)}}^{(1)},\boldsymbol{a}_{1}^{(2)},\cdots,\boldsymbol{a}_{J^{(2)}}^{(2)}$
as interventions decided on before the trial. By the definition of
$\boldsymbol{a}_{j}^{(2)}$, the two terms $\boldsymbol{U}_{1,n}^{(g)}$ and
${\boldsymbol{U}_{2\\_2,n}^{(g)\;*}}$ are independent. Thus, regular GEE
theory applies here, and Theorem 2 (Asymptotic Normality) follows. A detailed
proof of Theorem 2 can be found in Appendix B.2.

Appendix C presents the proofs for Theorem 1 and Theorem 2 using a log link
function as a concrete example.

##  5 Hypothesis testing

In a LAGO study, one of the main objectives is to evaluate the null hypothesis
of no overall intervention effect. This null hypothesis, represented by
$H_{0}:\boldsymbol{\beta}_{1}=0$, can be tested using standard methods, such
as the Chi-squared test. The validity of the Chi-squared test for this purpose
is guaranteed by Theorem 2.

In a c-LAGO study, a group indicator $R$ can be used to identify the
intervention group ($R=1$) and the control group ($R=0$). Let $\mu_{0}$ and
$\mu_{1}$ represent the mean outcome values in the control and intervention
groups, respectively. An alternative test for $H_{0}:\boldsymbol{\beta}_{1}=0$
is to test the implied $\tilde{H}_{0}:\mu_{0}=\mu_{1}$. Under the null
hypothesis, the stage 1 and stage 2 outcomes are independent because the stage
2 intervention has no effect on the stage 2 outcomes. Under the null
hypothesis, the distribution of the outcomes is the same regardless of the
intervention, and we can use any standard 1-degree-of-freedom test (such as
the Z-test) to compare the distribution in the intervention and the control
group. This leads to a test with the usual level of significance (either exact
or asymptotic, depending on the test used).

One can also take center characteristics $\boldsymbol{z}$ into account when
testing $H_{0}:\boldsymbol{\beta}_{1}=0$ in a c-LAGO study. The test involves
considering the null hypothesis that $\gamma=0$ in the model
$g\left(E\left(Y|\gamma,\boldsymbol{z};\boldsymbol{\beta}\right)\right)=\beta_{0}+\boldsymbol{\beta}_{2}^{T}\boldsymbol{z}+\gamma
R$. Under the null hypothesis, as before, there is independence between the
stages, and $\boldsymbol{\beta}_{1}=0$ implies that $\gamma=0$.

##  6 Confidence sets and confidence bands

The confidence set for the optimal intervention package is a list of
intervention package compositions that can be expected to include the optimal
intervention in 95$\%$ of such studies. To construct this confidence set, we
create a confidence interval for
$\mu=E\left(Y|\boldsymbol{x},\tilde{\boldsymbol{z}};\boldsymbol{\beta}^{*}\right)$
for a given value $\tilde{\boldsymbol{z}}$ and for each possible value of
$\boldsymbol{x}$. We first calculate a 95$\%$ confidence interval for
$g\left(\mu\right)$ as
$CI_{g(\mu)}=\left(1\;\boldsymbol{x}^{T}\;\tilde{\boldsymbol{z}}^{T}\right)\hat{\boldsymbol{\beta}}\pm
1.96\;\sqrt{\sigma^{2}_{\mu}(\hat{\boldsymbol{\beta}};\boldsymbol{x},\tilde{\boldsymbol{z}})}$
where
$\sigma^{2}_{\mu}(\hat{\boldsymbol{\beta}};\boldsymbol{x},\tilde{\boldsymbol{z}})$
can be calculated based on Theorem 2:
$\sigma^{2}_{\mu}(\hat{\boldsymbol{\beta}};\boldsymbol{x},\tilde{\boldsymbol{z}})=\left(1\;\boldsymbol{x}^{T}\;\tilde{\boldsymbol{z}}^{T}\right)n^{-1}\hat{J}\bigl{(}\hat{\boldsymbol{\beta}}\bigl{)}^{-1}\hat{V}\bigl{(}\hat{\boldsymbol{\beta}}\bigl{)}\hat{J}\bigl{(}\hat{\boldsymbol{\beta}}\bigl{)}^{-1}\left(1\;\boldsymbol{x}^{T}\;\tilde{\boldsymbol{z}}^{T}\right)^{T}$.
It follows that the 95$\%$ confidence interval for
$E\left(Y|\boldsymbol{x},\tilde{\boldsymbol{z}};\boldsymbol{\beta}^{*}\right)$
is $CI_{\mu}=g^{-1}\left(CI_{g(\mu)}\right)$. The confidence set can then be
calculated as
$CS\left(\boldsymbol{x}^{opt}\right)=\\{\boldsymbol{x}:CI_{\mu}\ni\theta\\}$.
That is, $CS\left(\boldsymbol{x}^{opt}\right)$ includes intervention packages
$\boldsymbol{x}$ for which $\theta$ (equation (2)) is inside the confidence
interval for the mean outcome under intervention package $\boldsymbol{x}$.
Because of Theorem 2, the confidence set $CS\left(\boldsymbol{x}^{opt}\right)$
contains $\boldsymbol{x}^{opt}$ with asymptotic probability at least 0.95.

Next, we construct confidence bands for the mean outcome under all different
intervention package compositions. These confidence bands have asymptotic
95$\%$ coverage simultaneously for all the intervention package compositions.
To calculate the confidence bands for the outcome
$E\left(Y|\boldsymbol{x},\tilde{\boldsymbol{z}};\boldsymbol{\beta}^{*}\right)$,
we first compute the 95$\%$ confidence bands for
$\left(1\;\boldsymbol{x}^{T}\;\tilde{\boldsymbol{z}}^{T}\right)\boldsymbol{\beta}^{*}$,
similar to Scheffe (1999) and Nevo, Lok and Spiegelman (2021) (Section 4 of
their Supplementary Material),

$CB_{g(\mu)}=\left(1\;\boldsymbol{x}^{T}\;\tilde{\boldsymbol{z}}^{T}\right)\hat{\boldsymbol{\beta}}\pm\sqrt{\chi_{0.95,p+q+1}^{2}\sigma^{2}_{\mu}(\hat{\boldsymbol{\beta}};\boldsymbol{x},\tilde{\boldsymbol{z}})},$
(8)

where $\chi_{0.95,p+q+1}^{2}$ is the 95th percentile of the $\chi_{p+q+1}^{2}$
distribution. $p+q+1$ is the dimension of $\boldsymbol{\beta}$, where $p$ and
$q$ are the dimensions of $\boldsymbol{a}$ and $\boldsymbol{z}$, respectively.
$CB_{\mu}=g^{-1}\bigl{(}CI_{g(\mu)}\bigl{)}$, and such confidence bands
guarantee asymptotic simultaneous 95$\%$ coverage for the mean outcome under
all intervention package compositions.

##  7 Simulations

Our methods were evaluated through simulation studies, which were organized
into four parts each based on 2000 simulated datasets. Simulation 1 considered
two scenarios of two-stage c-LAGO designs. Simulation 2 closely mimicked the
BetterBirth study as if a c-LAGO design was used. All non-adaptive design
parameters were taken directly from the BetterBirth study. Simulation 3
compared the performance of c-LAGO, uv-LAGO, and factorial designs. Simulation
4 evaluated the performance of three designs: c-LAGO, factorial design, and
MOST, employing the identity link function in the model for the continuous
outcome. The primary focus of Simulation 4 was on the power of each study
design.

In simulation 1 scenario 1, the data was generated with the same number of
centers $J$ in both the intervention and the control arm. The number of
centers in the intervention and the control arm was $J=6,10,20$, with
$n_{j}^{(1)}=50,100$, and $n_{j}^{(2)}=100,200$. The intervention consisted of
two components, $\boldsymbol{x}=(x_{1},x_{2})$, with the minimum and maximum
values of $x_{1}$ and $x_{2}$ being $\left[L_{1},U_{1}\right]=[0,2]$, and
$\left[L_{2},U_{2}\right]=[0,8]$, respectively. The recommended interventions
for stage 1 were set to be the middle of the range for the two intervention
components, $\boldsymbol{x}^{(1)}=(1,4)$. The model for the continuous outcome
was
$g\left(E\left(Y_{ij}|\boldsymbol{a}_{j},\boldsymbol{z}_{j}\;;\boldsymbol{\beta}\right)\right)=\boldsymbol{\beta}_{1}^{T}\boldsymbol{a}_{j}+\boldsymbol{\beta}_{2}^{T}\boldsymbol{z}_{j},$
where $g()$ is the logit link function. The logit link function was chosen as
this link function was used in the analysis of the BetterBirth study of
Section 8 to restrict the expected fractions to values between 0 and 1. The
true coefficient values for the intervention components
$\boldsymbol{\beta}_{1}^{*}$ were set as $(\beta_{11}^{*},\beta_{12}^{*})$
based on the estimated parameter values and their confidence intervals in the
final analysis of the BetterBirth study (Table 5): $(0.1863,0.15)$,
$(0.0438,0.17)$, $(0.1,0.2133)$, and $(0.1062,0.16)$. The exact true
coefficient values were chosen to facilitate the confidence set and confidence
bands calculations. More specifically,
$(\beta_{11}^{*}=0.1863,\beta_{12}^{*}=0.15)$ was chosen so that the true
optimal intervention was $(1,8)$ within three decimal places. Other
$\beta_{11}^{*}$ and $\beta_{12}^{*}$ values were chosen for the same reason.
A baseline center characteristic $Z\sim N(0,1)$ was also included with a true
coefficient value of $\beta_{z}^{*}=-0.2$. The aim was for a mean outcome
$\theta=0.8$, and the optimization problem was solved as described in equation
(2) to obtain the recommended interventions. For simplicity, no intercept was
included in the models.

In simulation 1 scenario 2, instead of having the same number of centers in
both stages, fewer centers and a smaller per-center sample size were included
in stage 1 than in stage 2. Specifically, the number of centers in stage 1 and
stage 2 were $J_{1}$=6, $J_{2}$=12 and $J_{1}$=10, $J_{2}$=20, respectively.
The per-center sample sizes were $n_{j}^{(1)}$ = 50 and $n_{j}^{(2)}$ = 200.
Other design parameters were the same as in simulation 1 scenario 1. We
incorporate variations in stage 1 interventions which accounts for the
possibility that centers may not adhere strictly to the interventions in stage
1. We do not incorporate variations from the recommended interventions in
stage 2, as those are based on previous stage outcomes and center specific.

Table 0a – 0c present selected results for simulation 1 scenario 1 and 2 using
a linear cost function with per unit cost for the two intervention components:
$C=(c_{1}=8,c_{2}=2)$. Complete results for simulation 1 with a linear cost
function can be found in Appendix E.1 \- E.3. Table 0a shows that for $J>6$,
both $\hat{{\beta}}_{11}$ and $\hat{{\beta}}_{12}$ had minimal relative bias.
The relative bias of $\hat{{\beta}}_{12}$ was smaller than that of
$\hat{{\beta}}_{11}$. The ratios between the mean of the estimated standard
error and the empirical standard error ranged from 0.7 – 0.9, however the
empirical coverage rates of the 95$\%$ confidence intervals for both
${{\beta}}_{11}$ and ${{\beta}}_{12}$ were close to 95$\%$. Despite having
thoroughly verified the simulation code’s accuracy, we cannot explain why
coverage remained satisfactory when the ratio between the mean of the
estimated standard error and the empirical standard error was small. However,
our analysis confirms that this phenomenon persisted.

Table 0b reports the bias and root mean squared errors for the estimated
optimal intervention components for a center with baseline center
characteristic ${z}$ equal to 0. The true optimal intervention components were
calculated using the method described in Remark 1 for linear cost functions,
with the true $\boldsymbol{\beta}$ values. The bias and root mean squared
error were small for both estimated optimal intervention components. The
estimated optimal intervention components based on stage 1 data (shown under
“Stage 1”) had higher bias and higher root mean squared error compared to the
final estimated optimal intervention based on all data (shown under “Stage
2/LAGO optimized”).

Table 0c displays information about the finite sample behavior of the
estimated optimal interventions, confidence sets, and confidence bands based
on stage 1 data and based on all data. The confidence set and the confidence
bands were computed by discretizing the two intervention components within
their lower and upper limits with the incrementation size of $0.1$. The steps
outlined in Section 6 were then applied. The true mean under the final
estimated intervention was close to 0.8 in most simulated datasets (see Table
0c column MenOpt2). The coverage rate for both the confidence set for the
optimal intervention package and the simultaneous confidence bands for the
intervention package components were very close to 95$\%$. The mean percentage
of the size of the confidence set as a percentage of the total sample space
ranged between 3$\%$ to 15$\%$ across different simulated datasets.

Table 1a – 1c present selected results for simulation 1 scenario 1 and 2 using
the cubic cost function. The cubic cost function for the two intervention
components $\boldsymbol{x}=(x_{1},x_{2})$ was defined as
$C(\boldsymbol{x})=0.05x_{1}^{3}-1.19x_{1}^{2}+10x_{1}+10+0.1x_{2}^{3}-0.7x_{2}^{2}+2x_{2}$.
The cubic cost function adapted the linear cost function from simulation 1 to
include an economy of scale at lower values of intervention component costs
and by including prohibitive cost as the intervention components neared their
upper limits. Complete results can be found in Appendix E.5. Table 1a
indicated minimal relative bias for both $\hat{{\beta}}_{11}$ and
$\hat{{\beta}}_{12}$. Similar to Table 0a, when the ratio between the mean of
the estimated standard error and the empirical standard error was small, the
empirical coverage rates remained satisfactory. Table 1b displays relatively
small bias and root mean squared error for the estimated optimal intervention
components based on all data (shown under “Stage 2/LAGO optimized”). The true
optimal intervention components were calculated using the method described in
Remark 1 for cubic cost functions, with the true $\boldsymbol{\beta}$ values.
The bias and root mean squared error were larger than those reported in Table
0b. Table 1c demonstrates satisfactory finite sample behavior in terms of the
estimated optimal intervention, confidence sets, and confidence bands. The
reported values aligned closely with those in Table 0c.

For simulation 2, we set all non-adaptive parameter values to those obtained
from the BetterBirth study, including the stage 1 interventions, the number of
centers, the per-center sample size, the intervention arm allocations, and the
distribution of the center-specific covariate, birth volume, $z$. For the
stage 2 and stage 3 interventions, we simulated a c-LAGO design. The results
of simulation 2 can be found in Appendix E.4.

In simulation 3, we compared the performance of the c-LAGO, uv-LAGO and
factorial designs. We set the number of centers $J=8$, and the per-center
sample sizes $n_{j}^{(1)}=n_{j}^{(2)}=10$ in each of the $K=2$ stages. The
intervention package consisted of two components,
$\boldsymbol{x}=(x_{1},x_{2})$ with minimum and maximum values of $x_{1}$ and
$x_{2}$ set to be the same as in simulation 1. The stage 1 outcomes were
simulated using a factorial design with interventions $(0,0)$, $(1,0)$,
$(0,4)$, and $(1,4)$, each with probability $1/4$. The stage 2 outcomes were
simulated using either a c-LAGO or a uv-LAGO design. The factorial design had
only one stage with the number of centers set to $J=16$. The model for the
continuous outcome was
$E\left(Y_{ij}|\boldsymbol{a}_{j},\boldsymbol{z}_{j};\boldsymbol{\beta}\right)=\beta_{0}+\boldsymbol{\beta}_{1}^{T}\boldsymbol{a}_{j}+\boldsymbol{\beta}_{2}^{T}\boldsymbol{z}_{j},$
and the errors were normally distributed with a standard deviation of 1.5. The
true parameter values were set to $\beta_{0}^{*}=0.1$, $\beta_{11}^{*}=0.2$,
$\beta_{12}^{*}=0.3$, and $\beta_{z}^{*}=-0.2$. The cost function for the two
intervention components was linear, with $C=(c_{1}=8,c_{2}=2)$ per unit. The
true optimal intervention for a center with $z=0$ was $x^{opt}=(0,3)$.

Under the null hypothesis of the c-LAGO design, with
$(\beta_{11}^{*},\beta_{12}^{*})=(0,0)$, the type-1 error of the component-
wise test (P degree of freedom Chi-squared test) was 4.7%, and 3.1% for the
two-sample means test. Under the null hypothesis of the uv-LAGO design, the
type-1 error of the component-wise test was 4.8%, and 3.6% for the two-sample
means test. Table 3 suggests that our approach was effective and led to
minimal finite sample bias, correct nominal coverage, and correct type-1
error.

In simulation 4, we evaluated the performance of the c-LAGO, factorial and
MOST designs with a focus on the power of the study design. The model for the
continuous outcome was a linear regression model with an identity link and
normal errors, with 3 stages and 100 participants per stage. Our goal was to
increase the expected outcome from 0.1 to 0.8 while minimizing cost, using a
linear cost function with $c_{1}=1$ and $c_{2}=5$. The outcome model was
$g(E(Y_{ij}|\boldsymbol{a}_{i};\boldsymbol{\beta}))=\beta_{0}+\beta_{11}a_{i1}+\beta_{12}a_{i2}$
with coefficients of $(\beta_{0},\beta_{11},\beta_{12})=(0.1,0.05,0.12)$ and
the errors were normally distributed with a standard deviation of 1.75. For
all designs, the stage 1 outcomes were generated according to a factorial
design with interventions (0,0), (2,0), (0,5), (2,5) each with a probability
of $1/4$. The factorial design continued with this intervention package
composition throughout. For the c-LAGO design, the stage 2 and stage 3
recommended interventions were determined with the goal of reaching an
expected mean outcome of 0.8 while minimizing cost and also achieving an
estimated power of 0.9. This power constraint was the main focus of our paper
under development and did not violate any of the assumptions outlined in
Section 2. In order to simplify and facilitate the comparison between c-LAGO,
factorial and MOST designs, the recommended interventions in c-LAGO were not
tailored to individual centers. The MOST design followed the same factorial
design for stage 1, then calculated the recommended intervention as in LAGO.
In the randomized controlled trial stage, half of the participants of MOST
were assigned to the control group and the other half to the treatment group.
The final analysis only included the outcomes after stage 1. The results in
Table 4 indicate that c-LAGO had minimal finite sample bias, correct nominal
coverage, and notably higher power compared to both the factorial design and
MOST.

Table 1: Selected simulation study results for simulation 1 scenario 1 and 2
with a linear cost function

$n_{j}^{(1)}$: number of patients in center $j$ at stage 1, $n_{j}^{(2)}$:
number of patients in center $j$ at stage 2.
$J$: number of centers for each stage.
%RelBias: percent relative bias
$100(\hat{\beta}-\beta^{\star})/\beta^{\star}$.
SE: mean estimated standard error, EMP.SD: empirical standard deviation.
CP95: empirical coverage rate of 95% confidence intervals.
Bias of $\hat{x}_{1}^{opt}$: bias of the first component of the estimated
optimal intervention, Bias of $\hat{x}_{2}^{opt}$: bias of the second
component of the estimated optimal intervention.
rMSE: root of mean squared
errors,$\Bigl{\\{}\operatorname{mean}\bigl{(}\bigl{\|}\hat{\boldsymbol{x}}^{opt}-\boldsymbol{x}^{opt}\bigl{\|}^{2}\bigl{)}\Bigl{\\}}^{1/2}$,
mean is taken over simulation iterations.
MeanOpt1: mean outcome under the stage 2 recommended intervention, calculated
using true coefficient values; MeanOpt2: mean outcome under the final
estimated optimal intervention based on all data, calculated using true
coefficient values.
$Q$2.5 and $Q$97.5: 2.5$\%$ and 97.5$\%$ quantiles.
SetCP95$\%$: empirical coverage percentage of confidence set for the optimal
intervention. SetPerc$\%$: mean percentage of the size of the confidence set
as a percent of the total sample space. BandsCP95$\%$: empirical coverage of
95$\%$ confidence band.

$\boldsymbol{\beta}^{*}=(\beta_{11}^{*},\beta_{12}^{*})$ | $n_{j}^{(1)}$ | $n_{j}^{(2)}$ | $J$ |  | $\hat{\boldsymbol{\beta}}_{11}$ |  |  | $\hat{\boldsymbol{\beta}}_{12}$ |
---|---|---|---|---|---|---|---|---|---
|  |  |  |  $\%$RelBias |  | $\frac{SE}{EMP.SD}$
---
($\times 100$)
CP95 |  $\%$RelBias |  | $\frac{SE}{EMP.SD}$
---
($\times 100$)
CP95
Scenario 1 ($J_{1}=J_{2}=J$) |  |  |  |  |  |
(0.1863, 0.15) | 50 | 100 | 6 | 1.64 | 84.8 | 96.1 | -0.46 | 76.1 | 95.0
|  |  | 10 | 1.78 | 90.6 | 94.9 | -0.55 | 79.3 | 94.9
|  |  | 20 | 1.54 | 97.6 | 95.1 | -0.37 | 91.2 | 95.2
|  | 200 | 6 | 1.44 | 77.9 | 95.2 | -0.42 | 70.9 | 95.5
|  |  | 10 | 1.91 | 86.1 | 95.4 | -0.49 | 73.6 | 95.5
|  |  | 20 | 1.14 | 97.8 | 95.4 | -0.26 | 88.2 | 94.8
Scenario 2a ($J_{1}=6,J_{2}=12$) |  |  |  |  |  |
(0.1863, 0.15) | 50 | 200 |  | 0.04 | 85.9 | 95.0 | -0.09 | 75.0 | 94.8
(0.1, 0.2133) | 50 | 200 |  | 0.56 | 93.7 | 95.8 | 0.14 | 84.3 | 95.7
(a)  Simulation study results for individual package component effects with a
linear cost function

$\boldsymbol{\beta}^{*}=(\tiny\beta_{11}^{*},\tiny\beta_{12}^{*})$ | $\boldsymbol{x}^{opt}$ | $n_{j}^{(1)}$ | $n_{j}^{(2)}$ |  | Stage 1 |  | Stage 2/LAGO optimized
---|---|---|---|---|---|---|---
|  |  |  |  |  | Bias of
---
$\hat{x}_{1}^{opt}$
($\times$100)
| Bias of
---
$\hat{x}_{2}^{opt}$
($\times$100)
| rMSE
---
($\times$100)
|  | Bias of
---
$\hat{x}_{1}^{opt}$
($\times$100)
| Bias of
---
$\hat{x}_{2}^{opt}$
($\times$100)
| rMSE
---
($\times$100)
Scenario 1 ($J_{1}=J_{2}=20$) |  |  |  |  |  |  |  |
(0.1863, 0.15) | (1,8) | 50 | 100 |  | 15.2 | 96.9 | 115.0 |  | -0.27 | 0.01 | 35.3
|  |  | 500 |  |  |  |  |  | -0.70 | 0.00 | 28.2
|  | 100 | 100 |  | 14.0 | 51.6 | 94.1 |  | 0.11 | 0.01 | 32.7
|  |  | 500 |  |  |  |  |  | -0.25 | 0.00 | 24.5
Scenario 2a ($J_{1}=6,J_{2}=12$) |  |  |  |  |  |  |  |
(0.1863, 0.15) | (1,8) | 50 | 200 |  | 20.7 | 341.5 | 191.1 |  | -1.5 | 1.5 | 52.3
(0.1, 0.2133) | (0, 6.5) | 50 | 200 |  | -46.2 | 239.2 | 172.9 |  | -0.1 | 0.1 | 39.2
(b) Simulation study results for estimated optimal intervention with a linear
cost function

$\boldsymbol{\beta}^{*}=(\tiny\beta_{11}^{*},\tiny\beta_{12}^{*})$ | $\boldsymbol{x}^{opt}$ | $n_{j}^{(1)}$ | $n_{j}^{(2)}$ |  | MeanOpt1
---
(Q2.5,Q97.5)
| MeanOpt2
---
(Q2.5,Q97.5)
| SetCP95
---
$\%$
| SetPerc
---
$\%$
| BandsCP95
---
$\%$
Scenario 1 ($J_{1}=J_{2}=20$) |  |  |  |  |
(0.1863, 0.15) | (1,8) | 50 | 100 | (0.655, 0.818) | (0.789, 0.811) | 95.3 | 5.5 | 95.2
|  |  | 500 |  | (0.794, 0.808) | 95.1 | 3.7 | 95.7
|  | 100 | 100 | (0.708, 0.816) | (0.791, 0.808) | 94.4 | 4.4 | 95.1
|  |  | 500 |  | (0.795, 0.805) | 95.2 | 2.8 | 95.8
Scenario 2a ($J_{1}=6,J_{2}=12$) |  |  |  |  |
(0.1863, 0.15) | (1,8) | 50 | 200 | (0.529, 0.826) | (0.773, 0.826) | 94.5 | 9.9 | 95.4
(0.1, 0.2133) | (0, 6.5) | 50 | 200 | (0.515, 0.857) | (0.783, 0.818) | 95.7 | 12.6 | 95.7
(c) Simulation study results for estimated optimal intervention, confidence
set and confidence band with a linear cost function

Table 2: Selected simulation study results for simulation 1 scenario 1 and 2
with a cubic cost function

$n_{j}^{(1)}$: number of patients in center $j$ at stage 1, $n_{j}^{(2)}$:
number of patients in center $j$ at stage 2.
$J$: number of centers for each stage.
%RelBias: percent relative bias
$100(\hat{\beta}-\beta^{\star})/\beta^{\star}$.
SE: mean estimated standard error, EMP.SD: empirical standard deviation.
CP95: empirical coverage rate of 95% confidence intervals.
Bias of $\hat{x}_{1}^{opt}$: bias of the first component of the estimated
optimal intervention, Bias of $\hat{x}_{2}^{opt}$: bias of the second
component of the estimated optimal intervention.
rMSE: root of mean squared
errors,$\Bigl{\\{}\operatorname{mean}\bigl{(}\bigl{\|}\hat{\boldsymbol{x}}^{opt}-\boldsymbol{x}^{opt}\bigl{\|}^{2}\bigl{)}\Bigl{\\}}^{1/2}$,
mean is taken over simulation iterations.
MeanOpt1: mean outcome under the stage 2 recommended intervention, calculated
using true coefficient values; MeanOpt2: mean outcome under the final
estimated optimal intervention based on all data, calculated using true
coefficient values.
$Q$2.5 and $Q$97.5: 2.5$\%$ and 97.5$\%$ quantiles.
SetCP95$\%$: empirical coverage percentage of confidence set for the optimal
intervention. SetPerc$\%$: mean percentage of the size of the confidence set
as a percent of the total sample space. BandsCP95$\%$: empirical coverage of
95$\%$ confidence band.

$\boldsymbol{\beta}^{*}=(\beta_{11}^{*},\beta_{12}^{*})$ | $n_{j}^{(1)}$ | $n_{j}^{(2)}$ | $J$ |  | $\hat{\boldsymbol{\beta}}_{11}$ |  |  | $\hat{\boldsymbol{\beta}}_{12}$ |
---|---|---|---|---|---|---|---|---|---
|  |  |  |  $\%$RelBias |  | $\frac{SE}{EMP.SD}$
---
($\times 100$)
CP95 |  $\%$RelBias |  | $\frac{SE}{EMP.SD}$
---
($\times 100$)
CP95
Scenario 1 ($J_{1}=J_{2}=J$) |  |  |  |  |  |
(0.1863, 0.15) | 50 | 100 | 6 | 2.31 | 81.9 | 95.8 | -0.79 | 78.1 | 95.5
|  |  | 10 | 4.59 | 95.8 | 95.4 | -1.63 | 88.3 | 95.0
|  |  | 20 | 4.35 | 101.3 | 95.7 | -1.49 | 94.1 | 96.0
|  | 200 | 6 | 1.36 | 83.9 | 95.5 | -0.52 | 78.8 | 95.5
|  |  | 10 | 2.68 | 94.2 | 95.9 | -0.89 | 85.6 | 96.0
|  |  | 20 | 3.90 | 96.3 | 96.4 | -1.29 | 88.1 | 96.1
Scenario 2a ($J_{1}=6,J_{2}=12$) |  |  |  |  |  |
(0.1863, 0.15) | 50 | 200 |  | 0.21 | 83.1 | 95.2 | -0.22 | 75.3 | 94.9
(0.1, 0.2133) | 50 | 200 |  | 2.54 | 88.1 | 95.6 | -0.31 | 74.5 | 95.8
(a) Simulation study results for individual package component effects with a
cubic cost function

$\boldsymbol{\beta}^{*}=(\beta_{11}^{*},\beta_{12}^{*})$ | $\boldsymbol{x}^{opt}$ | $n_{j}^{(1)}$ | $n_{j}^{(2)}$ |  | Stage 1 |  | Stage 2/LAGO optimized
---|---|---|---|---|---|---|---
|  |  |  |  |  | Bias of
---
$\hat{x}_{1}^{opt}$
($\times$100)
| Bias of
---
$\hat{x}_{2}^{opt}$
($\times$100)
| rMSE
---
($\times$100)
|  | Bias of
---
$\hat{x}_{1}^{opt}$
($\times$100)
| Bias of
---
$\hat{x}_{2}^{opt}$
($\times$100)
| rMSE
---
($\times$100)
Scenario 1 ($J_{1}=J_{2}=20$) |  |  |  |  |  |  |  |
(0.1863, 0.15) | (1.5,7.4) | 50 | 100 |  | 40.0 | 81.9 | 118.8 |  | 4.9 | 4.8 | 65.9
|  |  | 500 |  |  |  |  |  | 4.4 | -0.8 | 56.2
|  | 100 | 100 |  | 32.7 | 45.8 | 103.4 |  | 2.6 | 5.9 | 63.9
|  |  | 500 |  |  |  |  |  | 3.8 | 0.0 | 56.2
Scenario 2a ($J_{1}=6,J_{2}=12$) |  |  |  |  |  |  |  |
(0.1863, 0.15) | (1.5,7.4) | 50 | 200 |  | 64.4 | 314.0 | 187.6 |  | 9.7 | 19.9 | 79.8
(0.1, 0.2133) | (0, 6.5) | 50 | 200 |  | -64.2 | 261.3 | 170.0 |  | -13.8 | 7.3 | 54.1
(b) Simulation study results for estimated optimal intervention with a cubic
cost function

$\boldsymbol{\beta}^{*}=(\beta_{11}^{*},\beta_{12}^{*})$ | $\boldsymbol{x}^{opt}$ | $n_{j}^{(1)}$ | $n_{j}^{(2)}$ |  | MeanOpt1
---
(Q2.5,Q97.5)
| MeanOpt2
---
(Q2.5,Q97.5)
| SetCP95
---
$\%$
| SetPerc
---
$\%$
| BandsCP95
---
$\%$
Scenario 1 ($J_{1}=J_{2}=20$) |  |  |  |  |
(0.1863, 0.15) | (1.5,7.4) | 50 | 100 | (0.656, 0.812) | (0.769, 0.809) | 95.8 | 8.3 | 96.1
|  |  | 500 |  | (0.786, 0.809) | 95.3 | 5.3 | 95.6
|  | 100 | 100 | (0.708, 0.816) | (0.776, 0.807) | 95.9 | 7.0 | 96.3
|  |  | 500 |  | (0.789, 0.806) | 95.7 | 4.6 | 95.9
Scenario 2a ($J_{1}=6,J_{2}=12$) |  |  |  |  |
(0.1863, 0.15) | (1.5,7.4) | 50 | 200 | (0.509, 0.814) | (0.735, 0.821) | 94.8 | 11.5 | 95.8
(0.1, 0.2133) | (0, 6.5) | 50 | 200 | (0.515, 0.823) | (0.775, 0.820) | 95.9 | 14.7 | 95.8
(c) Simulation study results for estimated optimal intervention with a cubic
cost function

Table 3:  Comparing relative bias of $\hat{\boldsymbol{\beta}}$, its mean
estimated standard error over empirical standard deviation, coverage
probability and power among c-LAGO, uv-LAGO, and factorial design with a
linear cost function.

|  | c-LAGO |  | uv-LAGO |  | factorial
---|---|---|---|---|---|---
|  |  $\%$RelBias | $\frac{SE}{EMP.SD}$ |  |  $\%$RelBias | $\frac{SE}{EMP.SD}$ |  |  $\%$RelBias | $\frac{SE}{EMP.SD}$
Stage 1 | $\beta_{0}$ | -0.79 | 100.40 |  | -0.79 | 100.40 |  |  |
| $\beta_{11}$ | -4.55 | 97.99 |  | -4.55 | 97.99 |  |  |
| $\beta_{12}$ | 0.62 | 97.56 |  | 0.62 | 97.56 |  |  |
Stage 2 | $\beta_{0}$ | -25.93 | 96.52 |  | -26.23 | 95.64 |  | -2.26 | 99.93
| $\beta_{11}$ | -6.68 | 97.57 |  | -7.87 | 98.10 |  | 1.82 | 100.02
| $\beta_{12}$ | 4.67 | 98.16 |  | 4.79 | 98.46 |  | 0.04 | 102.46
|  | CP95
---
($\beta_{0}$,$\beta_{11}$,$\beta_{12}$)
(0.94, 0.94, 0.94) |  | (0.96, 0.98, 0.98) |  | (0.94, 0.95, 0.96)
| Power | (0.99, 0.93) |  | (1.00, 0.93) |  | (0.99, 0.91)
| (component-wise test, two-sample means test)

$\%$RelBias: percent relative bias
$100(\hat{\beta}-\beta^{\star})/\beta^{\star}$.
SE: mean estimated standard error, EMP.SD: empirical standard deviation.
CP95: empirical coverage rate of 95% confidence intervals.
component-wise test: P degree of freedom Chi-squared test.

Table 4:  Comparing bias of $\hat{\boldsymbol{\beta}}$, its mean estimated
standard error, coverage probability and power among c-LAGO, MOST and
factorial designs with a linear cost function.

|  | c-LAGO |  | factorial
---|---|---|---|---
|  | $\hat{\beta}_{0}$ | $\hat{\beta}_{11}$ | $\hat{\beta}_{12}$ |  | $\hat{\beta}_{0}$ | $\hat{\beta}_{11}$ | $\hat{\beta}_{12}$
Stage 1 | Bias | -0.004 | -0.004 | 0.002 | Bias | -0.004 | -0.001 | 0.001
| SE | 0.302 | 0.174 | 0.070 | SE | 0.175 | 0.101 | 0.040
Stage 1-2 | Bias | -0.011 | -0.027 | 0.005 | CP95 | 0.957 | 0.951 | 0.946
| SE | 0.185 | 0.114 | 0.051 | Power | 0.529 |  |
Stage 1-3 | Bias | -0.012 | -0.035 | 0.005 |  |  |  |
| SE | 0.145 | 0.094 | 0.041 |  |  |  |
| CP95 | 0.943 | 0.939 | 0.936 |  |  |  |
| Power | 0.871 |  |  |  |  |  |
|  | MOST |
| Optimization phase:
---
n=100
RCT: n=200
| Intervention
---
package effect:
Median (IQR)
| 0.555
---
(0.339, 0.730)
| C.P. intervention
---
package effect:
0.947 |
|  |  |  |  | Power | 0.489 |
| Optimization phase:
---
n=148*
RCT: n=152
| Intervention
---
package effect:
Median (IQR)
| 0.606
---
(0.437, 0.769)
| C.P. intervention
---
package effect:
0.949 |
|  |  |  |  |  | Power | 0.465 |

*: the sample size for the optimization phase must be divisible by 4 as it is a factorial design with 4 intervention groups.
SE: mean estimated standard error.
Intervention package effect: true mean under the intervention in the RCT.
C.P. intervention package effect: empirical coverage rate of 95% confidence
intervals for the treatment effect of the actual intervention in the RCT.
CP95: empirical coverage rate of 95% confidence intervals.
Power: two sample means test.
MOST Optimization trial n = 100: 1952 out of 2000 trials moved on to an RCT
(have at least 1 estimated coefficient greater than 0), 1092 out of 1952 RCTs
include only intervention component 1, 860 out of 1952 RCTs include only
intervention component 2, 0 out of 1952 RCTs include both intervention
components.
MOST Optimization trial n = 148: 1984 out of 2000 trials moved on to an RCT
(have at least 1 estimated coefficient greater than 0), 1165 out of 1952 RCTs
include only intervention component 1, 819 out of 1952 RCTs include only
intervention component 2, 0 out of 1952 RCTs include both intervention
components.

##  8 Illustrative example: The BetterBirth Study

We illustrated the LAGO design for continuous outcomes with the BetterBirth
study (Hirschhorn et al., 2015; Semrau et al., 2017). The BetterBirth study
aimed to improve maternal and neonatal health outcomes in Uttar Pradesh, India
by implementing the World Health Organization’s (WHO) Safe Childbirth
Checklist (SCC). This WHO checklist encouraged birth attendants to use
essential birth practices (EBPs) known to prevent complications at various
pause points during the delivery process.

The BetterBirth study consisted of three stages, the first two of which were
pilot studies. Stage 3 was a cluster randomized trial. The intervention
package used in stage 2 was modified based on feedback from stage 1, and
adjusted again in stage 3. Stage 1 included 2 centers, stage 2 included 4
centers, and stage 3 included a control and an intervention arm with 15
centers each. Outcome data was collected before and after the intervention was
implemented in the two pilot stages, and data was collected both before and
after the implementation of the intervention package for 5 centers in the
cluster randomized trial.

The outcome of interest was the proportion of EBPs performed out of all
possible birth practices measured during each stage, and we modeled this as a
continuous outcome. Births with 0 measured EBPs were excluded from the study
because whether EBPs were measured depends on the availability of researchers
to document them. The number of EBPs measured at each stage were 14, 19, and
18, respectively. The average proportion of EBPs performed out of all possible
birth practices measured for the three stages were 0.33, 0.28, and 0.42,
respectively.

We included the approximate monthly birth volume as a baseline center
characteristic $\boldsymbol{z}$. To avoid multicollinearity, we only
considered two out of the four intervention components: the duration of the
on-site intervention launch (in days) and the number of coaching visits after
the initial intervention launch, truncated to 40 visits or less. Because the
outcome of interest was the mean proportion of EBPs performed, we fit this
LAGO design by a GLM with a logit link function to restrict the expected
fractions to values between 0 and 1:

$logit\left(E\left(Y_{ij}|\boldsymbol{A}_{j}=\boldsymbol{a}_{j},\boldsymbol{X}_{j}=\boldsymbol{x}_{j},\boldsymbol{z}_{j};\boldsymbol{\beta}\right)\right)=\beta_{0}+\boldsymbol{\beta}_{1}^{T}\boldsymbol{a}_{j}+\boldsymbol{\beta}_{2}^{T}\boldsymbol{z}_{j}.$

Table 5 reports the estimated effects of the intervention package components
after each stage, based on all data available at the end of that stage. The
final analysis (last column) indicates that both the duration of the on-site
intervention launch and the number of coaching visits had positive effects.
The estimated effect of the number of coaching visits was also highly
significant across the different stages, and both the two sample means test
and the component-wise Chi-squared test had p-values less than 0.001.

Table 5: The BetterBirth Study: package component effect estimates, 95$\%$
confidence intervals after each stage, and estimated optimal intervention
package after each stage

|  | Stage 1
---
$n^{(1)}=113$
$\hat{\boldsymbol{\beta}}$ ($95\%$ CI)
| Stage 1-2
---
$n^{(1)}+n^{(2)}=2256$
$\hat{\boldsymbol{\beta}}$ ($95\%$ CI)
| Stage 1-3
---
$n^{(1)}+n^{(2)}+n^{(3)}$
$=7342$
$\hat{\boldsymbol{\beta}}$ ($95\%$ CI)
Intercept | 2.72 (1.20, 4.25) | -0.61 (-0.69, -0.53) | -0.138 (-0.156, -0.120)
Launch Duration (days) | -0.09 (-0.33, 0.14) | -0.003 (-0.11, 0.10) | 0.17 (0.11, 0.22)
Coaching Visits (per 5 visits) | 0.90 (0.79, 1.01) | 0.32 (0.31, 0.34) | 0.172 (0.167, 0.176)
Birth Volume (monthly, per 100) | -3.39 (-4.59, -2.18) | -0.166 (-0.180, -0.153) | -0.202 (-0.210, -0.195)
| $\hat{\boldsymbol{x}}^{opt}$
---
(using the linear cost function)
(1, 16) | (1, 36) | (5, 31)
| $\hat{\boldsymbol{x}}^{opt}$
---
(using the cubic cost function)
(1, 15.66) | (1, 35.38) | (4.73, 31.84)

CI: based on sandwich estimator for VAR($\hat{\boldsymbol{\beta}}$) (see
Theorem 2).
For (optimal) interventions, the first component is launch duration and the
second component is number of coaching visits.
The optimal intervention reported is for a center with an average birth volume
($z=175$).

Next, we found the optimal intervention package that results in a mean of
performed EBPs greater than 0.8 ($\theta=0.8$) while minimizing cost. Let
$x_{1}$ be the launch duration (in days) and $x_{2}$ be the number of coaching
visits, we used both the linear cost function:
$C_{1}(\boldsymbol{x})=800x_{1}+170x_{2}$ (Nevo, Lok and Spiegelman, 2021) and
the cubic cost function:
$C_{2}(\boldsymbol{x})=1700x_{1}-950x_{1}^{2}+220x_{1}^{3}+380x_{2}-24x_{2}^{2}+0.6x_{2}^{3}$
for the analysis. The cubic cost function revised the linear cost function to
include an economy of scale at lower values of intervention component costs
and by including prohibitive cost as the intervention components neared their
upper limits. Constraints were set such that $1\leq x_{1}\leq 5$ and $1\leq
x_{2}\leq 40$. For a center with an average birth volume ($z=175$), the
estimated optimal intervention package under the linear cost function
comprised a launch duration of 5 days and 31 coaching visits, with a total
cost of $\$9270$. With the cubic cost function, the estimated optimal
intervention package consisted of 3.97 days for launch duration and 35.50
coaching visits, at a total cost of $\$15629.04$. The closest integer values
for the estimated optimal intervention that can lead to a mean of performed
EBPs greater than 0.8 while minimizing cost are 4 days for launch duration and
36 coaching visits.

To determine the 95$\%$ confidence set $CS(\boldsymbol{x}^{opt})$ for the
optimal intervention, we examined a grid of all possible values of the
intervention components. The launch duration was incremented in steps of 0.01
days, ranging from 1 to 5 days (i.e., 1, 1.01, 1.02, …, 5), the number of
coaching visits was also incremented by 0.01, varying from 1 to 40 (i.e., 1,
1.01, 1.02, …, 40). The final 95$\%$ confidence set consisted of $1.01\times
10^{5}$ ($6.48\%$) out of the total $1.56\times 10^{6}$ intervention packages,
and included various combinations of the intervention package components.
Employing the linear cost function, the first, second, and third quartiles of
the cost within the confidence set were $\$8930$, $\$9150$, and $\$9380$,
respectively. Conversely, utilizing the cubic cost function, the corresponding
quartiles of the cost within the confidence set were $\$15356$, $\$16449$, and
$\$17481$. Under the linear cost function, with the estimated optimal
intervention as $\hat{\boldsymbol{x}}^{opt}=(5,31)$ (based on the confidence
bands for the proportion of EBPs performed under all possible intervention
package components), the 95% confidence interval for the proportion of EBPs
was (0.766, 0.834). Under the cubic cost function, with the estimated optimal
intervention as $\hat{\boldsymbol{x}}^{opt}=(3.97,35.50)$, the 95% confidence
interval for the proportion of EBPs was (0.780, 0.819).

##  9 Discussion

The LAGO design allows for changes in the composition of the intervention
package based on accumulating data from an ongoing trial at pre-specified
stages. LAGO could help prevent failed trials by adapting and optimizing the
intervention package composition while the trial is ongoing. LAGO is useful
for implementation trials, pragmatic trials, and clinical trials of
combination regimens. The methods described in this paper further increase the
flexibility of the LAGO design, and we anticipate that LAGO will be widely
adopted in intervention trials of combination interventions and combination
implementation strategies.

The objectives of LAGO studies are to determine the optimal intervention
package that achieves a pre-specified effect at minimal cost, test its
efficacy on the outcome of interest, and evaluate its effects. Variation in
the interventions is needed for LAGO to identify the treatment effect
parameters, which is an inherent feature of uv-LAGO and frequently observed in
large-scale public health studies. Furthermore, c-LAGO often recommends
different interventions for different centers based on their center
characteristics, thereby introducing variation. We have proven that the LAGO
design leads to consistent, asymptotically normal estimators when the outcomes
are continuous. The simulation studies show that LAGO has good finite sample
properties at reasonable sample sizes.

While LAGO designs in their current form do not encompass interim hypothesis
testing, it is indeed feasible to include futility stops in LAGO designs. By
defining a baseline acceptable power level, futility stops can be seamlessly
integrated into a LAGO trial. Specifically, if we fail to identify any
recommended intervention within the feasible bounds that is projected to yield
an adequate power after stage $k$, the trial may be terminated early. The type
I error does not increase from futility stops, since there is no strong
conclusion when a trial stops for futility (Snapinn et al., 2006).

As in Nevo, Lok and Spiegelman (2021), we used fixed baseline center
characteristics to account for any random center effects. Future research
could focus on developing LAGO for studies with center effects. In large-scale
intervention studies, cluster randomized trials are common, and the random
effects model is the standard approach for modeling these studies (Bell,
Fairbrother and Jones, 2019). Another potential extension is allowing centers
to participate in more than one stage of a LAGO study, which is likely to
occur in real world studies. In addition, future research could focus on
developing LAGO for individual-level interventions. Another important area for
future research is to systematically determine the optimal values of $K$,
$J_{k}$, $n_{jk}$, and the function $f$ of Remark 1, that provides the
recommended intervention in later stages.

The LAGO design will be applied to the PULESA-Uganda trial (Strengthening the
Blood Pressure Care and Treatment Cascade for Ugandans Living with HIV -
ImpLEmentation Strategies to SAve Lives; NIH UG3HL154501). Chronic HIV is a
risk factor for cardiovascular disease (CVD), and hypertension is the most
important driver of CVD risk. Blood pressure measurement is a key step in the
management of hypertension. PULESA will first explore current practice,
routines, barriers, and facilitators of evidence-based blood pressure care in
HIV clinical settings, and then use these findings to design an implementation
strategy to improve HIV-hypertension care. PULESA will determine the
effectiveness of the proposed intervention package and evaluate the
implementation strategy’s economic sustainability in Uganda’s Kampala and
Wakiso districts.

[Acknowledgments] The authors are extremely grateful to the participants of
the BetterBirth Study. The authors also thank Dr. Katherine Semrau for her
valuable contributions in explaining the BetterBirth Study data.

Judith J. Lok was supported by NSF Grant DMS 1854934. The BetterBirth Study
was supported by the Bill and Melinda Gates Foundation. The article contents
are the sole responsibility of the authors and do not necessarily represent
the official views of the Bill and Melinda Gates Foundation or the NSF.

## References

  * Beidas et al. (2022) [author] Beidas, Rinad SR. S., Dorsey, ShannonS., Lewis, Cara CC. C., Lyon, Aaron RA. R., Powell, Byron JB. J., Purtle, JonathanJ., Saldana, LisaL., Shelton, Rachel CR. C., Stirman, Shannon WiltseyS. W. and Lane-Fall, Meghan BM. B. (2022). Promises and pitfalls in implementation science from the perspective of US-based researchers: learning from a pre-mortem. Implementation Science 17 55\.
  * Bell, Fairbrother and Jones (2019) [author] Bell, AndrewA., Fairbrother, MalcolmM. and Jones, KelvynK. (2019). Fixed and random effects models: making an informed choice. Quality & Quantity 53 1051–1074.
  * Collins, Murphy and Strecher (2007) [author] Collins, Linda ML. M., Murphy, Susan AS. A. and Strecher, VictorV. (2007). The multiphase optimization strategy (MOST) and the sequential multiple assignment randomized trial (SMART): new methods for more potent eHealth interventions. American journal of preventive medicine 32 S112–S118.
  * Collins, Nahum-Shani and Almirall (2014) [author] Collins, Linda ML. M., Nahum-Shani, InbalI. and Almirall, DanielD. (2014). Optimization of behavioral dynamic treatment regimens based on the sequential, multiple assignment, randomized trial (SMART). Clinical Trials 11 426–434.
  * Collins et al. (2011) [author] Collins, Linda ML. M., Baker, Timothy BT. B., Mermelstein, Robin JR. J., Piper, Megan EM. E., Jorenby, Douglas ED. E., Smith, Stevens SS. S., Christiansen, Bruce AB. A., Schlam, Tanya RT. R., Cook, Jessica WJ. W. and Fiore, Michael CM. C. (2011). The multiphase optimization strategy for engineering effective tobacco use interventions. Annals of behavioral medicine 41 208–226.
  * Eisenberg and Gan (1983) [author] Eisenberg, BennettB. and Gan, Shi XinS. X. (1983). Uniform convergence of distribution functions. Proceedings of the American Mathematical Society 88 145–146.
  * Fahrmeir and Tutz (2013) [author] Fahrmeir, LudwigL. and Tutz, GerhardG. (2013). Multivariate statistical modelling based on generalized linear models. Springer Science & Business Media.
  * FDA (2016) [author] FDA (2016). Adaptive designs for medical device clinical studies: Guidance for industry and food and drug administration staff.
  * FDA (2019) [author] FDA (2019). Adaptive Design Clinical Trials for Drugs and Biologics Guidance for Industry. 2019.
  * Fogel (2018) [author] Fogel, David BD. B. (2018). Factors associated with clinical trials that fail and opportunities for improving the likelihood of success: a review. Contemporary clinical trials communications 11 156–164.
  * Greer (2010) [author] Greer, MonicaM. (2010). Electricity cost modeling calculations. Academic Press.
  * Hirschhorn et al. (2015) [author] Hirschhorn, Lisa RuthL. R., Semrau, KatherineK., Kodkany, BhalaB., Churchill, RobynR., Kapoor, AtulA., Spector, JonathanJ., Ringer, SteveS., Firestone, RebeccaR., Kumar, VishwajeetV. and Gawande, AtulA. (2015). Learning before leaping: integration of an adaptive study design process prior to initiation of BetterBirth, a large-scale randomized controlled trial in Uttar Pradesh, India. Implementation science 10 1–9.
  * Liang and Zeger (1986) [author] Liang, Kung-YeeK.-Y. and Zeger, Scott LS. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika 73 13–22.
  * McCullagh and Nelder (2019) [author] McCullagh, PeterP. and Nelder, John AJ. A. (2019). Generalized linear models. Routledge.
  * Nevo, Lok and Spiegelman (2021) [author] Nevo, DanielD., Lok, Judith JJ. J. and Spiegelman, DonnaD. (2021). Analysis of “Learn-As-you-GO”(LAGO) studies. The Annals of Statistics 49 793–819.
  * Scheffe (1999) [author] Scheffe, HenryH. (1999). The analysis of variance 72\. John Wiley & Sons.
  * Semrau et al. (2017) [author] Semrau, Katherine EAK. E., Hirschhorn, Lisa RL. R., Marx Delaney, MeganM., Singh, Vinay PV. P., Saurastri, RajivR., Sharma, NarenderN., Tuller, Danielle ED. E., Firestone, RebeccaR., Lipsitz, StuartS., Dhingra-Kumar, NeelamN. et al. (2017). Outcomes of a coaching-based WHO safe childbirth checklist program in India. New England Journal of Medicine 377 2313–2324.
  * Snapinn et al. (2006) [author] Snapinn, StevenS., Chen, Mon-GyM.-G., Jiang, QiQ. and Koutsoukos, TonyT. (2006). Assessment of futility in clinical trials. Pharmaceutical Statistics: The Journal of Applied Statistics in the Pharmaceutical Industry 5 273–281.
  * Stensland et al. (2014) [author] Stensland, Kristian DK. D., McBride, Russell BR. B., Latif, AsmaA., Wisnivesky, JuanJ., Hendricks, RyanR., Roper, NitinN., Boffetta, PaoloP., Hall, Simon JS. J., Oh, William KW. K. and Galsky, Matthew DM. D. (2014). Adult cancer clinical trials that fail to complete: an epidemic? JNCI: Journal of the National Cancer Institute 106\.
  * Van der Vaart (2000) [author] Van der Vaart, Aad WA. W. (2000). Asymptotic statistics 3\. Cambridge university press.

##  Appendix A  Overview of the Appendix

The Appendix is organized as follows: Appendix B presents the proofs of both
Theorem 1 and Theorem 2 for continuous outcome LAGO GLM with independent
errors and a general link function. Appendix C describes the proofs of both
Theorem 1 and Theorem 2 for LAGO GLM with independent errors and a log link
function. Appendix D extends the LAGO theory to include more than two stages.
Appendix E explains the motivation for the cubic cost function, and provides a
complete set of simulation results for the ones presented in Section 7 of the
main text.

##  Appendix B Proofs of Theorem 1 and Theorem 2 with general link function.

###  B.1 Proof of Theorem 1 for LAGO GLM with general link function:
consistency of $\hat{\boldsymbol{\beta}}$

First, the estimating equations from (3) of the main text are:

$\displaystyle 0=\boldsymbol{U}^{(g)}(\boldsymbol{\beta})$ (9)
$\displaystyle=\frac{1}{n}\bigg{[}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left(\frac{\partial}{\partial\boldsymbol{\beta}}E\left(Y_{ij}^{(1)}\mid\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}\right)\right)\left(Y_{ij}^{(1)}-E\left(Y_{ij}^{(1)}\mid\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}\right)\right)$
$\displaystyle\qquad+\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left(\frac{\partial}{\partial\boldsymbol{\beta}}E\left(Y_{ij}^{\left(2,n^{(1)}\right)}\mid\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}\right)\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\left(Y_{ij}^{\left(2,n^{(1)}\right)}-E\left(Y_{ij}^{\left(2,n^{(1)}\right)}\mid\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}\right)\right)\bigg{]}$
$\displaystyle=\frac{1}{n}\bigg{[}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)$
$\displaystyle\qquad+\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(Y_{ij}^{\left(2,n^{(1)}\right)}-g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\bigg{]}.$

Recall that $\boldsymbol{\beta}^{*}$ is the true value of $\boldsymbol{\beta}$
and let

$\displaystyle\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$ (10)
$\displaystyle=\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right]$
$\displaystyle+\sum_{j=1}^{J^{(2)}}\alpha_{j2}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})-g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right],$

where

$g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})=E\left(Y_{ij}^{\left(2\right)}\mid\boldsymbol{a}_{j}^{\left(2\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}\right).$

To show consistency of the estimator $\hat{\boldsymbol{\beta}}$, we show that
in spite of the fact that equation (9) does not consist of i.i.d. terms,
Theorem $5.9$ of Van der Vaart (2000) can be used. We show that its two
conditions are satisfied. First, we show

$\sup_{\boldsymbol{\beta}}\|\boldsymbol{U}^{(g)}(\boldsymbol{\beta})-\boldsymbol{u}^{(g)}(\boldsymbol{\beta})\|\xrightarrow{P}0.$
(11)

From equation (9) and equation (10), it follows that

$\displaystyle\boldsymbol{U}^{(g)}(\boldsymbol{\beta})-\boldsymbol{u}^{(g)}(\boldsymbol{\beta})$
(12)
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right]$
$\displaystyle+\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(Y_{ij}^{\left(2,n^{(1)}\right)}-g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right.\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.\left.\pm\;g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right)\right]$
$\displaystyle-\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right]$
$\displaystyle-\sum_{j=1}^{J^{(2)}}\alpha_{j2}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})-g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right]$
$\displaystyle=G_{1,n}^{(g)}+G_{2,n}^{(g)}+G_{3,n}^{(g)}+G_{4,n}^{(g)}+G_{5,n}^{(g)},$

where

$\displaystyle G_{1,n}^{(g)}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})\right)\right],$
(13)

$\displaystyle
G_{2,n}^{(g)}=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(2,n^{(1)})}-g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};\boldsymbol{\beta}^{*})\right)\right],$
(14)

$\displaystyle G_{3,n}^{(g)}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};\boldsymbol{\beta}^{*})\right)\right.$
(15)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\left.-\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right)\right],$

$\displaystyle G_{4,n}^{(g)}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right.$
(16)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\left.-\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right],$

$\displaystyle G_{5,n}^{(g)}=$
$\displaystyle\sum_{j=1}^{J^{(1)}}\left[\left(\alpha_{j1}-\frac{n_{j}^{(1)}}{n}\right)\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right]$
(17)
$\displaystyle+\sum_{j=1}^{J^{(1)}}\left[\left(\frac{n_{j}^{(1)}}{n}-\alpha_{j1}\right)\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})\right)\right]$
$\displaystyle+\sum_{j=1}^{J^{(2)}}\left[\left(\alpha_{j2}-\frac{n_{j}^{(2)}}{n}\right)\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right]$
$\displaystyle+\sum_{j=1}^{J^{(2)}}\left[\left(\frac{n_{j}^{(2)}}{n}-\alpha_{j2}\right)\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right)\right].$

The five terms $G_{1,n}^{(g)}$, $G_{2,n}^{(g)}$, $G_{3,n}^{(g)}$,
$G_{4,n}^{(g)}$, and $G_{5,n}^{(g)}$ will be discussed separately. We show the
supremum over $\boldsymbol{\beta}$ of each term converges to 0 in probability,
then the triangle inequality implies equation (11).

For $G_{1,n}^{(g)}$ from equation (13), we show that
$\sup_{\boldsymbol{\beta}}\left\|G_{1,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$ by using the concept of Donsker classes from empirical
process theory. Let
$O^{(1)}_{ij}=(Y^{(1)}_{ij},\boldsymbol{a}^{(1)}_{j},\boldsymbol{z}^{(1)}_{j})$
be the observed data for patient $i$ from center $j$ in stage 1, and let

$\Psi_{\boldsymbol{\beta}}(O^{(1)}_{ij})=\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})\right).$

We show that the class of functions
$\mathcal{F}=\left(\Psi_{\boldsymbol{\beta}}:\boldsymbol{\beta}\in\boldsymbol{\mathcal{B}}\right)$
is Donsker by showing that

$\left\|\Psi_{\boldsymbol{\beta}_{1}}(O^{(1)}_{ij})-\Psi_{\boldsymbol{\beta}_{2}}(O^{(1)}_{ij})\right\|\leq
m(O^{(1)}_{ij})\left\|\boldsymbol{\beta}_{1}-\boldsymbol{\beta}_{2}\right\|\quad\text{
for every }\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2},$ (18)

for some measurable function $m$ with $E\left(m^{2}\right)<\infty$; that
$\mathcal{F}$ is Donsker then follows from Example 19.7 in Van der Vaart
(2000). By applying the Mean Value Theorem to each row of
$\Psi_{\boldsymbol{\beta}_{1}}(O^{(1)}_{ij})-\Psi_{\boldsymbol{\beta}_{2}}(O^{(1)}_{ij})$
separately, we find that

$\displaystyle\Psi_{\boldsymbol{\beta}_{1}}(O^{(1)}_{ij})-\Psi_{\boldsymbol{\beta}_{2}}(O^{(1)}_{ij})$
$\displaystyle=\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\Psi_{\boldsymbol{\beta}}(O^{(1)}_{ij})\right)\left(\boldsymbol{\beta}_{1}-\boldsymbol{\beta}_{2}\right),$

where for each row of
$\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\Psi_{\boldsymbol{\beta}}(O^{(1)}_{ij})$,
$\tilde{\boldsymbol{\beta}}$ may take a different value between
$\boldsymbol{\beta}_{1}$ and $\boldsymbol{\beta}_{2}$. By Assumption 1 and 4
from the main text, $g()$ and $\Psi_{\boldsymbol{\beta}}(O^{(1)}_{ij})$ are
continuously differentiable and both $\boldsymbol{\beta}$ and $O^{(1)}_{ij}$
take values in a compact space, each element in
$\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\Psi_{\boldsymbol{\beta}}(O^{(1)}_{ij})$
is bounded. It follows that

$\displaystyle\left\|\Psi_{\boldsymbol{\beta}_{1}}(O^{(1)}_{ij})-\Psi_{\boldsymbol{\beta}_{2}}(O^{(1)}_{ij})\right\|$
$\displaystyle=\left\|\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\Psi_{\boldsymbol{\beta}}(O^{(1)}_{ij})\right)(\boldsymbol{\beta}_{1}-\boldsymbol{\beta}_{2})\right\|$
$\displaystyle\leq
C_{1}^{(g)}\left\|\boldsymbol{\beta}_{1}-\boldsymbol{\beta}_{2}\right\|,$

where $C_{1}^{(g)}$ is a constant. Equation (18) follows, so the class of
functions
$\mathcal{F}=\left(\Psi_{\boldsymbol{\beta}}:\boldsymbol{\beta}\in\boldsymbol{\mathcal{B}}\right)$
is Donsker. By Theorem 19.4 (Van der Vaart (2000)), $\mathcal{F}$ is also
Glivenko-Cantelli. Since for each $\boldsymbol{\beta}$,
$E\left(\Psi_{\boldsymbol{\beta}}\left(O_{ij}^{(1)}\right)\right)=0$, then

$\sup_{\boldsymbol{\beta}}\left\|\frac{1}{n_{j}^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\Psi_{\boldsymbol{\beta}}\left(O_{ij}^{(1)}\right)-E\left(\Psi_{\boldsymbol{\beta}}\left(O_{ij}^{(1)}\right)\right)\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0.$

Notice that the $O_{ij}^{(1)}$ are i.i.d. for each $j$ separately,
${n_{j}^{(1)}}/{n}<1$, and $j$ is finite, so

$\sup_{\boldsymbol{\beta}}\left\|\frac{n_{j}^{(1)}}{n}\frac{1}{n_{j}^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\Psi_{\boldsymbol{\beta}}\left(O_{ij}^{(1)}\right)\right\|\leq\sup_{\boldsymbol{\beta}}\left\|\frac{1}{n_{j}^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\Psi_{\boldsymbol{\beta}}\left(O_{ij}^{(1)}\right)\right\|.$

We conclude that
$\sup_{\boldsymbol{\beta}}\left\|G_{1,n}^{(g)}\right\|\xrightarrow{P}0.$

For the term $G_{2,n}^{(g)}$ from equation (14), let $Y_{ij}^{(2)}$ be the
(counterfactual) outcomes under $\boldsymbol{a}_{j}^{(2)}$ and let
$\epsilon_{ij}^{(2)}$ be the corresponding errors that patient $i$ in center
$j$ would have experienced under intervention $\boldsymbol{a}_{j}^{(2)}$. We
derive

$\displaystyle G_{2,n}^{(g)}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(2,n^{(1)})}-g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};\boldsymbol{\beta}^{*})\right)\right]$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\epsilon_{ij}^{(2,n^{(1)})}\right]$

By Assumption 5 from the main text, replacing $\epsilon_{ij}^{(2,n^{(1)})}$ by
the new error terms $\epsilon_{ij}^{(2)}$ does not change the distribution of
$G_{2,n}^{(g)}$, so it suffices to show that for

$\displaystyle\widetilde{G}_{2,n}^{(g)}=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\epsilon_{ij}^{(2)}\right],$

$\sup_{\boldsymbol{\beta}}\left\|\widetilde{G}_{2,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$.

$\displaystyle\widetilde{G}_{2,n}^{(g)}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\pm\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right]\epsilon_{ij}^{(2)}$
$\displaystyle=\widetilde{G}_{2\\_1,n}^{(g)}+\widetilde{G}_{2\\_2,n}^{(g)},$

where

$\widetilde{G}_{2\\_1,n}^{(g)}=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\epsilon_{ij}^{(2)}\right],$
$\widetilde{G}_{2\\_2,n}^{(g)}=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)-\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right]\epsilon_{ij}^{(2)}.$
(19)

We show that both
$\sup_{\boldsymbol{\beta}}\left\|\widetilde{G}_{2\\_1,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$ and
$\sup_{\boldsymbol{\beta}}\left\|\widetilde{G}_{2\\_2,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$. By the triangle inequality, then
$\sup_{\boldsymbol{\beta}}\left\|\widetilde{G}_{2,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$.

Let
$O^{(2)}_{ij}=(Y^{(2)}_{ij},\boldsymbol{a}^{(2)}_{j},\boldsymbol{z}^{(2)}_{j})$
be the counterfactual data for patient $i$ from center $j$ in stage 2 under
$\boldsymbol{a}^{(2)}_{j}$ and $\boldsymbol{z}^{(2)}_{j}$ and let

$\Psi_{\boldsymbol{\beta}}^{(2)}(O^{(2)}_{ij})=\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\epsilon_{ij}^{(2)}.$

Following the same argument as for $G_{1,n}^{(g)}$, for fixed value of $j$,
the class of functions
$\mathcal{F}_{2}=\left(\Psi_{\boldsymbol{\beta}}^{(2)}:\boldsymbol{\beta}\in\boldsymbol{\mathcal{B}}\right)$
is a Donsker class and
$\sup_{\boldsymbol{\beta}}\|\widetilde{G}_{2\\_1,n}^{(g)}\|\xrightarrow{P}0$.

For each of the $P$ components of $\widetilde{G}_{2\\_2,n}^{(g)}$ from
equation (19), by the Mean Value Theorem, we derive

$\displaystyle\widetilde{G}_{2\\_2,\;p}^{(g)}=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left[\left.\frac{\partial}{\partial
a_{p}}\right|_{\tilde{\boldsymbol{a}}_{jp}(\boldsymbol{\beta})}\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{a},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right]\left(\boldsymbol{A}_{jp}^{(2,n^{(1)})}-\boldsymbol{a}_{jp}^{(2)}\right)\epsilon_{ij}^{(2)}\right\\}.$

By Assumption 4 from the main text and Lemma B.1 from the Appendix, it follows
that

$\displaystyle\sup_{\boldsymbol{\beta}}\left\|\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)-\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right]\epsilon_{ij}^{(2)}\right\|$
$\displaystyle\leq\sup_{\boldsymbol{\beta}}\max_{j}{\left\|\left.\frac{\partial}{\partial
a_{p}}\right|_{\tilde{\boldsymbol{a}}_{jp}(\boldsymbol{\beta})}\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{a},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right\|}\max_{j}\left\|\left(\boldsymbol{A}_{jp}^{(2,n^{(1)})}-\boldsymbol{a}_{jp}^{(2)}\right)\right\|\sup_{ij}{\left\|\epsilon_{ij}^{(2)}\right\|}$
$\displaystyle\xrightarrow{P}0.$

This argument can be applied to all $P$ components of
$\widetilde{G}_{2\\_2,n}^{(g)}$, so
$\sup_{\boldsymbol{\beta}}\|\widetilde{G}_{2\\_2,n}^{(g)}\|\xrightarrow{P}0$.
Hence
$\sup_{\boldsymbol{\beta}}\left\|G_{2,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$.

Consider $G_{3,n}^{(g)}$ from equation (15). By the Mean Value Theorem, the
supremum over $\boldsymbol{\beta}$ of each of the $P$ components of
$G_{3,n}^{(g)}$ is

$\displaystyle\sup_{\boldsymbol{\beta}}\left\|\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};\boldsymbol{\beta}^{*})\right)\right.\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\left.\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right)\right]\rule{0.0pt}{22.76228pt}\right\|$
$\displaystyle=\sup_{\boldsymbol{\beta}}\left\|\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left.\frac{\partial}{\partial
a_{p}}\right|_{\tilde{\boldsymbol{a}}_{jp}(\boldsymbol{\beta})}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{a},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right]\right.\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\Biggl{.}\Biggl{.}\left(\boldsymbol{A}_{jp}^{(2,n^{(1)})}-\boldsymbol{a}_{jp}^{(2)}\right)\Biggr{\\}}\Biggr{\|}$
$\displaystyle\leq\sup_{\boldsymbol{\beta}}\max_{j}{\left\|\left.\frac{\partial}{\partial
a_{p}}\right|_{\tilde{\boldsymbol{a}}_{jp}(\boldsymbol{\beta})}\left(\frac{\partial}{\partial\boldsymbol{\beta}_{p}}g^{-1}(\boldsymbol{a},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right\|}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\max_{j}\biggr{\|}\left(\boldsymbol{A}_{jp}^{(2,n^{(1)})}-\boldsymbol{a}_{jp}^{(2)}\right)\biggr{\|}$
$\displaystyle\xrightarrow{P}0.$

The convergence to 0 follows from Assumption 4 of the main text and Lemma B.1
of the Appendix, so
$\sup_{\boldsymbol{\beta}}\left\|G_{3,n}^{(g)}\right\|\xrightarrow{P}0$. By
applying the same argument to $G_{4,n}^{(g)}$ (equation (16)), we conclude
that also
$\sup_{\boldsymbol{\beta}}\left\|G_{4,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$.

Consider the first term of $G_{5,n}^{(g)}$ (equation (17)), by Assumption 4
from the main text and $\alpha_{jk}=\lim_{n\rightarrow\infty}n_{j}^{(k)}/n$,
for the finite many $j$, it follows that

$\displaystyle\sup_{\boldsymbol{\beta}}\left\|\sum_{j=1}^{J^{(1)}}\left[\left(\alpha_{j1}-\frac{n_{j}^{(1)}}{n}\right)\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right]\right\|$
$\displaystyle\leq$
$\displaystyle\sup_{\boldsymbol{\beta}}\max_{j}\left\|\left(\alpha_{j1}-\frac{n_{j}^{(1)}}{n}\right)\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right\|$
$\displaystyle\xrightarrow{P}$ $\displaystyle\;0.$

By applying the same argument to the other terms of $G_{5,n}^{(g)}$ and the
triangle inequality,
$\sup_{\boldsymbol{\beta}}\left\|G_{5,n}^{(g)}\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$.

Thus,

$\displaystyle\sup_{\boldsymbol{\beta}}\left\|\boldsymbol{U}^{(g)}(\boldsymbol{\beta})-\boldsymbol{U}^{(g)}(\boldsymbol{\beta})\right\|$
$\displaystyle\leq\sup_{\boldsymbol{\beta}}\left\|G_{1,n}^{(g)}\right\|+\sup_{\boldsymbol{\beta}}\left\|G_{2,n}^{(g)}\right\|+\sup_{\boldsymbol{\beta}}\left\|G_{3,n}^{(g)}\right\|+\sup_{\boldsymbol{\beta}}\left\|G_{4,n}^{(g)}\right\|+\sup_{\boldsymbol{\beta}}\left\|G_{5,n}^{(g)}\right\|$
$\displaystyle\xrightarrow{P}0.$

The second condition in Theorem 5.9 of Van der Vaart (2000) is

$\inf_{\boldsymbol{\beta}:\left\|\boldsymbol{\beta}-\boldsymbol{\beta}^{\star}\right\|>0}\|\boldsymbol{u}^{(g)}(\boldsymbol{\beta})\|>0=\left\|\boldsymbol{u}^{(g)}\left(\boldsymbol{\beta}^{\star}\right)\right\|.$

Equation (10) implies that
$\left\|\boldsymbol{u}^{(g)}\left(\boldsymbol{\beta}^{\star}\right)\right\|=0$.
Furthermore, $\boldsymbol{u}^{(g)}\left(\boldsymbol{\beta}\right)$ is the same
as for a fixed two-stage design with
$a_{1}^{(1)},\cdots,a_{J^{(1)}}^{(1)},a_{1}^{(2)},\cdots,a_{J^{(2)}}^{(2)}$ as
interventions decided on before the trial, so regular GEE theory applies here.
In order for LAGO to work properly, we need variations in the intervention
components to identify the treatment effect parameter. The uniqueness of
$\boldsymbol{\beta}^{*}$ as a maximizer or zero has been studied by various
authors, see e.g. Chapter 2.2 of Fahrmeir and Tutz (2013). Thus, provided
there is enough variation in the intervention, the second condition in Theorem
5.9 of Van der Vaart (2000) is often also satisfied and we conclude that
$\hat{\boldsymbol{\beta}}$ is consistent.

###  B.2 Proof of Theorem 2 for LAGO GLM with general link function:
asymptotic normality of $\hat{\boldsymbol{\beta}}$

To prove asymptotic normality of the final estimator
$\hat{\boldsymbol{\beta}}$, we first prove

######  Lemma B.1.

Under Assumption 3, there exist $\boldsymbol{a}_{j}^{(2)}$, which is the
probability limit of $\boldsymbol{A}_{j}^{(2,n^{(1)})}$ as
$n^{(1)}\rightarrow\infty$ such that for stage $2$,
$\max_{j=1,\cdots,J^{(2)}}\left\|\boldsymbol{A}_{j}^{(2,n^{(1)})}-\boldsymbol{a}_{j}^{(2)}\right\|\xrightarrow{P}0$.

###### Proof.

By Assumption 3 and the Continuous Mapping Theorem, for each $j$,
$\boldsymbol{A}_{j}^{(2,n^{(1)})}\xrightarrow{P}\boldsymbol{a}_{j}^{(2)}$.
Since the number of centers in stage $2$ is fixed, taking the maximum over all
$j$ is a continuous operation. Lemma B.1 follows by the Continuous Mapping
Theorem. ∎

Under Assumptions 1 $-$ 5, and Lemma B.1, we show that

$\footnotesize{\sqrt{n}\left(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}\right)\xrightarrow{D}N(0,J\left(\boldsymbol{\beta^{*}}\right)^{-1}V\left(\boldsymbol{\beta^{*}}\right)J\left(\boldsymbol{\beta^{*}}\right)^{-1}),}$
(20)

where the explicit forms of $J\left(\boldsymbol{\beta^{*}}\right)$ and
$V\left(\boldsymbol{\beta^{*}}\right)$ are given in Theorem 2. The
corresponding estimators are

$\displaystyle\hat{J}(\hat{\boldsymbol{\beta}})$
$\displaystyle=\sum_{j=1}^{J^{(1)}}\frac{n_{j}^{(1)}}{n}\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\hat{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}+\sum_{j=1}^{J^{(2)}}\frac{n_{j}^{(2)}}{n}\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\hat{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2},$ (21) $\displaystyle\hat{V}(\hat{\boldsymbol{\beta}})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\left\\{\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\hat{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\sum_{i=1}^{n_{j}^{(1)}}\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\hat{\boldsymbol{\beta}})\right)^{2}\right\\}$
(22)
$\displaystyle+\frac{1}{n}\sum_{j=1}^{J^{(2)}}\left\\{\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\hat{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{A}_{j}^{(2,n^{(1)})},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2}\sum_{i=1}^{n_{j}^{(2)}}\left(Y_{ij}^{(2,n^{(1)})}-g^{-1}(\boldsymbol{A}_{j}^{(2,n^{(1)})},\boldsymbol{z}_{j}^{(2)};\hat{\boldsymbol{\beta}})\right)^{2}\right\\}.$

Applying the Mean Value Theorem to each component of
$\boldsymbol{U}^{(g)}(\boldsymbol{\beta})$ from equation (3),

$\footnotesize{0=\boldsymbol{U}^{(g)}(\hat{\boldsymbol{\beta}})=\boldsymbol{U}^{(g)}(\boldsymbol{\beta}^{*})+\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}({\boldsymbol{\beta}})\right)(\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}^{*}})^{T},}$

where for each row of
$\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}({\boldsymbol{\beta}})$,
$\tilde{\boldsymbol{\beta}}$ takes a possibly row-dependent value on the line
between $\hat{\boldsymbol{\beta}}$ and $\boldsymbol{\beta}^{*}$. It follows
that

$\footnotesize{\sqrt{n}\left(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}\right)=-\sqrt{n}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}\left({\boldsymbol{\beta}}\right)\right)^{-1}\boldsymbol{U}^{(g)}\left(\boldsymbol{\beta}^{*}\right).}$
(23)

We first show that
$\left(-\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}\left({\boldsymbol{\beta}}\right)\right)$
converges in probability to $J\left(\boldsymbol{\beta^{*}}\right)$. Then we
show that $\sqrt{n}\;\boldsymbol{U}^{(g)}\left(\boldsymbol{\beta}^{*}\right)$
converges to a normal distribution with mean 0 and variance
$V\left(\boldsymbol{\beta^{*}}\right)$. Equation (20) then follows from
Slutsky’s Theorem.

Because of equation (3) of the main text,

$\displaystyle\frac{\partial}{\partial\boldsymbol{\beta}}\boldsymbol{U}^{(g)}({\boldsymbol{\beta}})=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\frac{\partial}{\partial\boldsymbol{\beta}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\right]$
(24)
$\displaystyle\qquad\qquad+\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\frac{\partial}{\partial\boldsymbol{\beta}}\left[\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\left(Y_{ij}^{\left(2,n^{(1)}\right)}-g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\right]$
$\displaystyle=G_{6,n}^{(g)}(\boldsymbol{\beta})+G_{7,n}^{(g)}(\boldsymbol{\beta}),$

where

$\displaystyle G_{6,n}^{(g)}(\boldsymbol{\beta})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left\\{\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)Y_{ij}^{(1)}\right.$
$\displaystyle\left.-\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})-\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\},$ $\displaystyle G_{7,n}^{(g)}(\boldsymbol{\beta})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)Y_{ij}^{(2,n^{(1)})}\right.$
$\displaystyle\left.-\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}.$

Let

$\displaystyle G_{6,n}^{(g)*}(\boldsymbol{\beta})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left\\{\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})\right.$
$\displaystyle\left.-\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})-\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\},$ $\displaystyle G_{7,n}^{(g)*}(\boldsymbol{\beta})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right.$
$\displaystyle\left.-\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})-\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}.$

We show that
$\frac{\partial}{\partial\boldsymbol{\beta}}\boldsymbol{U}^{(g)}(\tilde{\boldsymbol{\beta}})-\left(G_{6,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)+G_{7,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)\right)$
converges in probability to 0. First,

$\footnotesize{\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}(\tilde{\boldsymbol{{\beta}}})-\left(G_{6,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)+G_{7,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)\right)=\left(G_{6,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{6,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)\right)+\left(G_{7,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)\right).}$

We show that both terms converge in probability to 0. For the first term,
using
$Y_{ij}^{(1)}=g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};{\boldsymbol{\beta}^{*}})+\epsilon_{ij}^{(1)}$,

$\displaystyle
G_{6,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{6,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)$
(25) $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left\\{\left(\left.\frac{\partial^{2}}{\partial{\boldsymbol{\beta}}^{2}}\right|_{\tilde{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};{\boldsymbol{\beta}})-\left.\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}\right|_{{\boldsymbol{\beta}^{*}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta}^{*})\right.$
$\displaystyle\left.+\left(\left.\frac{\partial^{2}}{\partial{\boldsymbol{\beta}^{2}}}\right|_{\tilde{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};{\boldsymbol{\beta}})\right)\;\epsilon_{ij}^{(1)}\right.+\left(\left.\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}\right|_{{\boldsymbol{\beta}^{*}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta}^{*})$
$\displaystyle-\left(\left.\frac{\partial^{2}}{\partial{\boldsymbol{\beta}^{2}}}\right|_{\tilde{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};{\boldsymbol{\beta}})\right)g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\tilde{\boldsymbol{\beta}})$
$\displaystyle+\left.\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{{\boldsymbol{\beta}^{*}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta})\right)^{\otimes
2}-\left(\left.\frac{\partial}{\partial{\boldsymbol{\beta}}}\right|_{\tilde{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};{\boldsymbol{\beta}})\right)^{\otimes
2}\right\\}.$

Notice that the $pq$-th entry of the first term of equation (25) equals

$\displaystyle\sum_{j=1}^{J^{(1)}}\frac{n_{j}^{(1)}}{n}\left\\{\left(\left.\frac{\partial}{\partial{\boldsymbol{\beta}}_{p}}\right|_{\tilde{\boldsymbol{\beta}}_{p}}\left.\frac{\partial}{\partial{\boldsymbol{\beta}}_{q}}\right|_{\tilde{\boldsymbol{\beta}}_{q}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};{\boldsymbol{\beta}})\right)g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta}^{*})\right.$
(26)
$\displaystyle\qquad\qquad\qquad-\left.\left(\left.\frac{\partial}{\partial{\boldsymbol{\beta}}_{p}}\right|_{{\boldsymbol{\beta}}^{*}_{p}}\left.\frac{\partial}{\partial{\boldsymbol{\beta}}_{q}}\right|_{{\boldsymbol{\beta}}^{*}_{q}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};\boldsymbol{\beta}^{*})\right\\}.$

Since $\hat{\boldsymbol{\beta}}$ and therefore $\tilde{\boldsymbol{\beta}}$ is
consistent, $\tilde{\boldsymbol{\beta}}\xrightarrow{P}\boldsymbol{\beta}^{*}$.
By Assumption 4 and the Continuous Mapping Theorem, equation (26) converges in
probability to 0 for each entry. Since $J^{(1)}$ is fixed, $g()$ is in
$C^{2}$, and ${n_{j}^{(1)}}/{n}\leq 1$, we conclude that equation (25) goes to
0 in probability. Next, following the same arguments as for ${G_{1,n}^{(g)}}$
in the proof of consistency (Appendix B.1),

$\footnotesize{\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left(\left.\frac{\partial^{2}}{\partial{\boldsymbol{\beta}}^{2}}\right|_{\tilde{\boldsymbol{\beta}}}g^{-1}(\boldsymbol{a}_{j}^{\left(1\right)},\boldsymbol{z}_{j}^{\left(1\right)};{\boldsymbol{\beta}})\right)\;\epsilon_{ij}^{(1)}\xrightarrow{P}0.}$

We conclude that indeed
$G_{6,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{6,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$.

For
$G_{7,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)$,
we first show
$\sup_{\boldsymbol{\beta}}\left\|G_{7,n}^{(g)}({\boldsymbol{\beta}})-G_{7,n}^{(g)*}\left(\boldsymbol{\beta}\right)\right\|\xrightarrow{P}0$,
after which we show that
$\left(G_{7,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)*}\left(\boldsymbol{\beta}^{*}\right)\right)\xrightarrow{P}0.$
Using that
$Y_{ij}^{(2,n^{(1)})}=g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};{\boldsymbol{\beta}^{*}})+\epsilon_{ij}^{(2,n^{(1)})}$,
from equation (24)

$\displaystyle
G_{7,n}^{(g)}(\boldsymbol{\beta})-G_{7,n}^{(g)^{*}}(\boldsymbol{\beta})=G_{7\\_1,n}^{(g)}(\boldsymbol{\beta})+G_{7\\_2,n}^{(g)}(\boldsymbol{\beta})+G_{7\\_3,n}^{(g)}(\boldsymbol{\beta})+G_{7\\_4,n}^{(g)}(\boldsymbol{\beta})$

where

$\footnotesize{G_{7\\_1,n}^{(g)}(\boldsymbol{\beta})=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)\epsilon_{ij}^{(2,n^{(1)})}\right\\},}$
$\displaystyle G_{7\\_2,n}^{(g)}(\boldsymbol{\beta})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};\boldsymbol{\beta}^{*})\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\left.-\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta}^{*})\right\\},$
$\displaystyle G_{7\\_3,n}^{(g)}(\boldsymbol{\beta})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\left.-\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}^{2}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right\\},$
$\displaystyle G_{7\\_4,n}^{(g)}(\boldsymbol{\beta})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left\\{\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2}-\left(\frac{\partial}{\partial\boldsymbol{\beta}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}.$

Following the same arguments as for $G_{2,n}^{(g)}$, $G_{3,n}^{(g)}$, and
$G_{4,n}^{(g)}$ in the proof of consistency (Appendix B.1), respectively, we
conclude that suprema over $\boldsymbol{\beta}$ of
$G_{7\\_1,n}^{(g)}(\boldsymbol{\beta})$,
$G_{7\\_2,n}^{(g)}(\boldsymbol{\beta})$,
$G_{7\\_3,n}^{(g)}(\boldsymbol{\beta})$ and
$G_{7\\_4,n}^{(g)}(\boldsymbol{\beta})$ converge to 0 in probability. Next,
notice that

$\displaystyle
G_{7,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)^{*}}(\boldsymbol{\beta}^{*})$
$\displaystyle=\left(G_{7,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)^{*}}(\tilde{\boldsymbol{\beta}})\right)+\left(G_{7,n}^{(g)^{*}}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)^{*}}\left(\boldsymbol{\beta}^{*}\right)\right).$

Since
$\sup_{\boldsymbol{\beta}}\left\|G_{7,n}^{(g)}(\boldsymbol{\beta})-G_{7,n}^{(g)^{*}}(\boldsymbol{\beta})\right\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$, it follows that
$G_{7,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)^{*}}(\tilde{\boldsymbol{\beta}})\stackrel{{\scriptstyle
P}}{{\rightarrow}}0.$ Similar to equation (25), we apply the Mean Value
Theorem to conclude that
$G_{7,n}^{(g)^{*}}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)^{*}}\left(\boldsymbol{\beta}^{*}\right)\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$. Combining, we conclude
$G_{7,n}^{(g)}(\tilde{\boldsymbol{\beta}})-G_{7,n}^{(g)^{*}}\left(\boldsymbol{\beta}^{*}\right)\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$.

It follows that
$\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}({\boldsymbol{{\beta}}})-\left(G_{6,n}^{(g)^{*}}(\boldsymbol{\beta}^{*})+G_{7,n}^{(g)^{*}}(\boldsymbol{\beta}^{*})\right)\xrightarrow{P}0$.
Thus, calculating the limit of
$-\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}({\boldsymbol{{\beta}}})$
is equivalent to calculating the limit of
$-\left(G_{6,n}^{(g)^{*}}(\boldsymbol{\beta}^{*})+G_{7,n}^{(g)^{*}}(\boldsymbol{\beta}^{*})\right)$.
We therefore show that
$-\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}({\boldsymbol{{\beta}}})$
converges in probability to $J(\boldsymbol{\beta}^{*})$ from Theorem 2.

$\displaystyle
G_{6,n}^{(g)^{*}}(\boldsymbol{\beta}^{*})-\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left\\{-\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left\\{-\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}+\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left\\{\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}$
$\displaystyle=\sum_{j=1}^{J^{(1)}}\frac{n_{j}^{(1)}}{n}\left\\{-\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}-\sum_{j=1}^{J^{(1)}}\alpha_{j1}\left\\{-\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}$ $\displaystyle\xrightarrow{P}0.$

Similarly,

$\footnotesize{G_{7,n}^{(g)^{*}}(\boldsymbol{\beta}^{*})-\sum_{j=1}^{J^{(1)}}\alpha_{j2}\left\\{-\left.\left(\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(2)},\boldsymbol{z}_{j}^{(2)};\boldsymbol{\beta})\right)^{\otimes
2}\right\\}\xrightarrow{P}0.}$

We conclude that

$\footnotesize{\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\tilde{\boldsymbol{\beta}}}\boldsymbol{U}^{(g)}({\boldsymbol{{\beta}}})\xrightarrow{P}J(\boldsymbol{\beta}^{*}).}$
(27)

To show that $\sqrt{n}\>\>\boldsymbol{U}\left(\boldsymbol{\beta}^{*}\right)$
from equation (23) converges in distribution to
$N\left(0,V\left(\boldsymbol{\beta}^{*}\right)\right)$, we derive

$\displaystyle\sqrt{n}\ \boldsymbol{U}^{(g)}(\boldsymbol{\beta}^{*})$ (28)
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=1}^{J^{(1)}}\sum_{i=1}^{n_{j}^{(1)}}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta})\right)\left(Y_{ij}^{(1)}-g^{-1}(\boldsymbol{a}_{j}^{(1)},\boldsymbol{z}_{j}^{(1)};\boldsymbol{\beta}^{*})\right)$
$\displaystyle\qquad+\frac{1}{\sqrt{n}}\sum_{j=1}^{J^{(2)}}\sum_{i=1}^{n_{j}^{(2)}}\left(\left.\frac{\partial}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}^{*}}g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};\boldsymbol{\beta})\right)\left(Y_{ij}^{\left(2,n^{(1)}\right)}-g^{-1}(\boldsymbol{A}_{j}^{\left(2,n^{(1)}\right)},\boldsymbol{z}_{j}^{\left(2\right)};\boldsymbol{\beta}^{*})\right)$