# Adaptive Randomization Methods for Sequential Multiple Assignment Randomized
Trials (SMARTs) via Thompson Sampling

Peter Norwood1,∗<EMAIL_ADDRESS>Marie Davidian2∗∗
<EMAIL_ADDRESS>and Eric Laber3∗∗∗
1Quantum Leap Healthcare Collaborative<EMAIL_ADDRESS>San Francisco
California  U.S.A.
2Department of Statistics  North Carolina State University  Raleigh  North
Carolina  U.S.A.
3Department of Statistical Science  Duke University  Durham  North Carolina
U.S.A

###### Abstract

Response-adaptive randomization (RAR) has been studied extensively in
conventional, single-stage clinical trials, where it has been shown to yield
ethical and statistical benefits, especially in trials with many treatment
arms. However, RAR and its potential benefits are understudied in sequential
multiple assignment randomized trials (SMARTs), which are the gold-standard
trial design for evaluation of multi-stage treatment regimes. We propose a
suite of RAR algorithms for SMARTs based on Thompson Sampling (TS), a widely
used RAR method in single-stage trials in which treatment randomization
probabilities are aligned with the estimated probability that the treatment is
optimal. We focus on two common objectives in SMARTs: (i) comparison of the
regimes embedded in the trial, and (ii) estimation of an optimal embedded
regime. We develop valid post-study inferential procedures for treatment
regimes under the proposed algorithms. This is nontrivial, as (even in single-
stage settings) RAR can lead to nonnormal limiting distributions of
estimators. Our algorithms are the first for RAR in multi-stage trials that
account for nonregularity in the estimand. Empirical studies based on real-
world SMARTs show that TS can improve in-trial subject outcomes without
sacrificing efficiency for post-trial comparisons.

###### keywords:

Inverse probability weighted estimator; Precision medicine; Response-adaptive
randomization, Treatment regime

##  1 Introduction

A treatment regime, also known as an adaptive treatment strategy, is a
sequence of decision rules, one for each key decision point in a patient’s
disease progression, that maps accumulated information on a patient to a
recommended treatment. (Tsiatis et al., 2020). Sequential multiple assignment
randomized trials (SMARTs) are the gold standard for the study of treatment
regimes (Lavori and Dawson, 2004; Murphy, 2005) and have been applied
successfully in a range of areas, including cancer, addiction, HIV/prevention,
and education (e.g., Kidwell, 2014; Bigirumurame et al., 2022; Lorenzoni et
al., 2023). A SMART involves multiple stages of randomization, each stage
corresponding to a decision point, where the sets of possible treatments may
depend on baseline and interim information.

Figure 1 depicts a two-stage SMART to evaluate behavioral intervention
strategies for cancer pain management (Somers et al., 2023). Subjects were
randomized at baseline with equal probability to one of two first-stage
interventions: Pain Coping Skills Training (PCST) with five sessions (PCST-
Full), or one session (PCST-Brief). At the end of the first stage, subjects
were classified as responders if they experienced a 30% reduction in their
pain score from baseline and as nonresponders otherwise. Subjects were then
assigned with equal probability to one of two second-stage interventions
depending on their first-stage intervention and response status. This design,
like any fixed (nonadaptive) SMART design, can be represented as a single-
stage trial in which subjects are randomized at baseline among a set of fixed
regimes known as the SMART’s embedded regimes (Tsiatis et al., 2020, Chapter
9). The cancer pain SMART has eight embedded regimes determined by the stage 1
treatment, stage 2 treatment for responders, and stage 2 treatment for
nonresponders, e.g., one embedded regime assigns PCST-Full initially followed
by no further treatment if the subject responds and another round of PCST-Full
otherwise. A key goal was to evaluate the embedded regimes on the basis of
mean percent reduction in pain from baseline at the end of stage 2 and to
identify a regime that yields the greatest reduction.

Figure 1:  Schematic depicting the design of the SMART for evaluation of
behavioral interventions for cancer pain management. The eight embedded
regimes implied by the design are of the form “Give $a$ initially; if
response, give $b$, otherwise if nonresponse give $c$,” where regimes 1,…, 8
correspond to $(a,b,c)=$ (0,0,1), (0,0,2), (0,1,2), (0,1,1), (1,3,4), (1,3,2),
(1,4,2), (1,4,4), respectively.

As in the preceding example, most SMARTs use fixed randomization probabilities
at each stage. Although fixed and balanced randomization in SMARTs yields high
power for comparing treatments and treatment regimes (Murphy, 2005), ethical
and pragmatic considerations suggest that updating randomization probabilities
based on accruing information can improve outcomes for trial subjects,
increase enrollment, and decrease dropout (Food and Drug Administration,
2019). Response-adaptive randomization (RAR) uses accumulating information to
skew randomization probabilities toward promising treatments and has long been
used in single-stage randomized clinical trials (RCTs, Kim et al., 2011;
Berry, 2015; Wang and Yee, 2019). There is a lengthy literature on RAR
methods; reviews for single-stage RCTs include Hu and Rosenberger (2006),
Berry et al. (2010), and Atkinson and Biswas (2013). RAR procedures are
typically formalized as a multi-arm bandit in which each treatment option is a
bandit “arm” (Berry and Fristedt, 1985; Villar et al., 2015; Lattimore and
Szepesvári, 2020); strong theoretical and empirical support for bandits has
grown from their use in applications such as mHealth and telehealth (Tewari
and Murphy, 2017; Liu et al., 2023). In single-stage RCTs, the advantages of
RAR are most pronounced in trials with a large number of treatments because
RAR algorithms oversample more favorable treatments and undersample less
favorable ones (Berry, 2011). RAR also presents challenges, e.g., early
evidence can lead the randomization to become “stuck” favoring suboptimal
treatments (Thall, Fox, and Wathen, 2015), and there have been spirited
debates about the potential pitfalls of RAR (Viele et al., 2020; Proschan and
Evans, 2020; Villar et al., 2021).

As the number of treatments and regimes evaluated in a SMART can be large,
integration of RAR into SMARTs could have significant benefits, yet
development is limited. Cheung, Chakraborty, and Davidson (2014) propose
SMART-AR, an RAR method based on Q-learning (e.g., Tsiatis et al., 2020,
Section 5.7.1) that adapts randomization probabilities so that treatments with
large estimated Q-functions are more likely to be selected. However, this
approach does not account for uncertainty in the estimated Q-functions nor the
potential information gain associated with each treatment, so the resulting
design may be inefficient. Wang, Wu, and Wahed (2022) propose RA-SMART, an RAR
method for two-stage SMARTs with the same treatments at each stage. The method
does not incorporate delayed effects, i.e., does not account for treatment
effects at stage 2 that depend on stage 1 treatment. Because the possibility
of interactions among treatments at different stages is key to the development
of treatment regimes, RAR schemes that acknowledge delayed effects are
desirable.

In this article, we propose RAR approaches for SMARTs based on Thompson
sampling (TS, Thompson, 1933), a popular bandit algorithm also known as
probability matching (Agrawal and Goyal, 2012; Russo et al., 2018). The idea
underlying TS is that a treatment’s randomization probability should be based
on confidence that the treatment is optimal. Classical TS measures this
confidence with a posterior probability; aligned with standard SMART
methodology, we adopt a frequentist perspective in which confidence is
assessed using a confidence distribution (Xie and Singh, 2013). The proposed
methods are applicable to SMARTs used to evaluate a fixed set of treatment
regimes and/or identify an optimal regime and, by construction, account for
delayed effects.

In Section 2, we present the statistical framework and a variant of TS for
SMARTs, and we use the latter in Section 3 to construct the proposed RAR
approaches. One class of methods randomizes subjects “up-front” to entire
embedded regimes, while the other sequentially randomizes subjects at each
decision point. We propose estimators for the marginal mean outcome under a
regime in Section 3 and argue in Section 4 that they are consistent and
asymptotically normal; this is nontrivial, as adaptive randomization can lead
to nonnormal limits for plug-in estimators (Hadad et al., 2021; Zhang et al.,
2021b). Simulations demonstrating performance are reported in Section 5.

##  2 Background and preliminaries

###  2.1 Notation and assumptions

We first review SMARTs with fixed (i.e., nonadaptive) randomization. Consider
a SMART with $K$ stages, where, at stage $k=1,\ldots,K$, $\mathcal{A}_{k}$ is
the (finite) set of available treatment options. Let $A_{k}\in\mathcal{A}_{k}$
be the treatment assigned at stage $k$, and let
$\overline{A}_{k}=(A_{1},\ldots,A_{k})\in\mathcal{A}_{1}\times\cdots\times\mathcal{A}_{k}$
denote all treatments given through stage $k$. Let $\boldsymbol{X}_{1}$ be a
set of baseline subject variables collected prior to administration of stage 1
treatment, and let $\boldsymbol{X}_{k}$ comprise variables collected between
stages $k-1$ and $k$, $k=2,\ldots,K$. Define
$\overline{\boldsymbol{X}}_{k}=(\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{k})$,
$k=1,\ldots,K$; and let $\boldsymbol{H}_{1}=\boldsymbol{X}_{1}$ and
$\boldsymbol{H}_{k}=(\boldsymbol{X}_{1},A_{1},\boldsymbol{X}_{2},A_{2},\dots,A_{k-1},\boldsymbol{X}_{k})=(\overline{\boldsymbol{X}}_{k},\overline{\boldsymbol{A}}_{k-1})$,
$k=2,\ldots,K$, denote the information at the time $A_{k}$ is assigned, with
$\mathcal{H}_{k}$ denoting the domain of $\boldsymbol{H}_{k}$, $k=1,\ldots,K$.
Let $Y$ be the real-valued outcome of interest, which is measured or
constructed from $\overline{\boldsymbol{X}}_{K}$ and information measured
after stage $K$ at some specified follow-up time, coded so that larger values
are more favorable. Elements of $\boldsymbol{X}_{k}$ may include current and
past measures of the patient’s health status, previous measures of the
outcome, and response status. For a subject with history
$\boldsymbol{H}_{k}=\boldsymbol{h}_{k}$ at stage $k$,
$\Psi_{k}(\boldsymbol{h}_{k})\subseteq\mathcal{A}_{k}$ is the set of feasible
treatment options in $\mathcal{A}_{k}$, where $\Psi_{k}$ maps
$\mathcal{H}_{k}$ to subsets of $\mathcal{A}_{k}$ (Tsiatis et al., 2020,
Section 6.2.2). E.g., at stage $k=2$ of the cancer pain SMART,
$\Psi_{2}(\boldsymbol{h}_{2})=\\{1,2\\}$ for a subject who does not respond to
treatment 0 given at the first stage. At stage $k$, a subject is randomized to
option $a_{k}\in\Psi_{k}(\boldsymbol{h}_{k})$ with probability
$P(A_{k}=a_{k}|\boldsymbol{H}_{k}=\boldsymbol{h}_{k})$; often,
$P(A_{k}=a_{k}|\boldsymbol{H}_{k}=\boldsymbol{h}_{k})=1/|\Psi_{k}(\boldsymbol{h}_{k})|$
for all $k$.

A decision rule $d_{k}:\mathcal{H}_{k}\rightarrow\mathcal{A}_{k}$ maps an
individual’s history to a recommended treatment at stage $k$, where
$d_{k}(\boldsymbol{h}_{k})\in\Psi_{k}(\boldsymbol{h}_{k})$ for all
$\boldsymbol{h}_{k}\in\mathcal{H}_{k}$. A treatment regime
$\mathbf{d}=\left(d_{1},\dots,d_{K}\right)$ is a sequence of decision rules,
where $d_{k}$ is the rule for stage $k$, $k=1,\ldots,K$. The mean outcome that
would be achieved if the population were to receive treatment according to
$\mathbf{d}$ can be characterized in terms of potential outcomes. For any
$\overline{\boldsymbol{a}}_{k-1}=(a_{1},\ldots,a_{k-1})\in\mathcal{A}_{k-1}$,
let $\boldsymbol{X}_{k}^{*}(\overline{\boldsymbol{a}}_{k-1})$ denote the
potential subject information that would accrue between stages $k-1$ and $k$,
and define
$\overline{\boldsymbol{X}}_{k}^{*}(\overline{\boldsymbol{a}}_{k-1})=\\{\boldsymbol{X}_{1},\boldsymbol{X}_{2}^{*}(a_{1}),\ldots,\boldsymbol{X}_{k}^{*}(\overline{\boldsymbol{a}}_{k-1})\\}$,
$k=2,\ldots,K$. The potential history at stage $k\geq 2$ is thus
$\boldsymbol{H}_{k}^{*}(\overline{\boldsymbol{a}}_{k-1})=\\{\overline{\boldsymbol{X}}_{k}^{*}(\overline{\boldsymbol{a}}_{k-1}),\overline{\boldsymbol{a}}_{k-1}\\}$,
with
$\boldsymbol{X}_{1}^{*}(\overline{\boldsymbol{a}}_{0})=\boldsymbol{H}_{1}^{*}(\overline{\boldsymbol{a}}_{0})=\boldsymbol{H}_{1}$.
For any $\overline{\boldsymbol{a}}_{K}\in\overline{\mathcal{A}}_{K}$, let
$Y^{*}(\overline{\boldsymbol{a}}_{K})$ be the potential outcome that would be
achieved under $\overline{\boldsymbol{a}}_{K}$. For $k=2,\ldots,K$, the
potential accrued information between stages $k-1$ and $k$ under regime
$\mathbf{d}$ is then
$\boldsymbol{X}^{*}_{k}(\mathbf{d})=\sum_{\overline{\boldsymbol{a}}_{k-1}\in\overline{\mathcal{A}}_{k-1}}\boldsymbol{X}^{*}_{k}(\overline{\boldsymbol{a}}_{k-1})\prod^{k-1}_{v=1}I[d_{v}\\{\boldsymbol{H}_{v}^{*}(\overline{\boldsymbol{a}}_{v-1})\\}=a_{v}]$,
and the potential outcome under $\mathbf{d}$ is
$Y^{*}(\mathbf{d})=\sum_{\overline{\boldsymbol{a}}_{K}\in\overline{\mathcal{A}}_{K}}Y^{*}(\overline{\boldsymbol{a}}_{K})\prod^{K}_{k=1}I[d_{k}\\{\boldsymbol{H}_{k}^{*}(\overline{\boldsymbol{a}}_{k-1})\\}=a_{k}]$,
where $I(\cdot)$ is the indicator function. The mean outcome for regime
$\mathbf{d}$, known as the value of $\mathbf{d}$, is
$\mathcal{V}(\mathbf{d})=E\\{Y^{*}(\mathbf{d})\\}$.

Define
$\mathcal{W}^{*}=\\{\boldsymbol{X}_{2}^{*}(a_{1}),\ldots,\boldsymbol{X}_{K}^{*}(\overline{\boldsymbol{a}}_{K-1}),Y^{*}(\overline{\boldsymbol{a}}_{K})$
for all $\overline{\boldsymbol{a}}_{K}\in\mathcal{A}_{K}\\}$ to be the set of
all potential outcomes. For a given regime $\mathbf{d}$, identification of
$\mathcal{V}(\mathbf{d})$ is possible from the data
$(\overline{\boldsymbol{X}}_{K},\overline{\boldsymbol{A}}_{K},Y)$ under the
following assumptions, which are discussed extensively elsewhere (e.g.,
Tsiatis et al., 2020, Section 6.2.4) and which we adopt: (i) consistency,
$Y=Y^{*}(\overline{\boldsymbol{A}}_{K})$,
$\boldsymbol{H}_{k}=\boldsymbol{H}_{k}^{*}(\overline{\boldsymbol{A}}_{k-1})$,
$k=1,\ldots,K$; (ii) positivity,
$P(A_{k}=a_{k}|\boldsymbol{H}_{k}=\boldsymbol{h}_{k})>0$ for all
$a_{k}\in\Psi_{k}(\boldsymbol{h}_{k})$ and all
$\boldsymbol{h}_{k}\in\mathcal{H}_{k}$, $k=1,\ldots,K$; and (iii) sequential
ignorability, $\mathcal{W}^{*}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}A_{k}\,|\,\boldsymbol{H}_{k}$, $k=1,\ldots,K$,
where “$\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}$” denotes statistical independence. We also
assume that there is no interference among individuals nor multiple versions
of a treatment. In a SMART with nonadaptive randomization, (ii) is true by
design, and (iii) is guaranteed by randomization; with RAR, (ii) holds if
constraints are imposed on the randomization probabilities (discussed shortly)
and (iii) holds because randomization probabilities are known features of the
history.

Denote the $m$ regimes embedded in a SMART as
$\mathbf{d}^{1},\ldots,\mathbf{d}^{m}$. A common primary analysis is the
comparison of $\mathcal{V}(\mathbf{d}^{j})$, $j=1,\ldots,m$, i.e., the mean
outcomes that would be achieved if the population were to receive treatments
according to $\mathbf{d}^{1},\ldots,\mathbf{d}^{m}$. Another common analysis
is identification of an optimal embedded regime, $\mathbf{d}^{\mathrm{opt}}$,
satisfying
$\mathcal{V}(\mathbf{d}^{\mathrm{opt}})\geq\mathcal{V}(\mathbf{d}^{j})$,
$j=1,\ldots,m$.

###  2.2 Subject accrual and progression processes

Ordinarily, it is assumed that subjects enroll in a clinical trial by a
completely random process over a planned accrual period. Under such a process,
subjects in a SMART progress through the stages in a staggered fashion. Thus,
at any point during the trial, enrolled subjects will have reached different
stages, with only some having the outcome $Y$ ascertained. By reaching stage
$k$, we mean that a subject has completed the previous stages and that $A_{k}$
has been assigned. E.g., in a SMART with $K=2$ stages, at any time, there may
be subjects who have received only a first-stage treatment, for whom
$(\boldsymbol{X}_{1},A_{1})$ is available; subjects who have reached stage 2
but have not completed follow up, for whom
$(\boldsymbol{X}_{1},A_{1},\boldsymbol{X}_{2},A_{2})$ is available; and
subjects who have completed the trial, for whom
$(\boldsymbol{X}_{1},A_{1},\boldsymbol{X}_{2},A_{2},Y)$ is available.

Ideally, randomization probabilities would be updated each time a new
randomization is needed. However, for logistical or other reasons, it may be
more feasible in practice to update the probabilities on a pre-set schedule.
For definiteness, we consider a SMART that will enroll $N$ subjects over $T$
weeks, with randomization probabilities updated weekly. We develop two
randomization schemes: the first, which we term up-front, randomizes subjects
to an embedded regime at baseline; the second, which we term sequential,
randomizes subjects at each stage. While these two schemes are equivalent in a
SMART with fixed randomization probabilities, the sequential scheme under RAR
allows a subject’s randomization probabilities to depend on information
collected as they progress through the trial.

We assume that a group of subjects of random size enrolls at each week $t$
and, depending on the randomization scheme, up-front or sequential, requires
assignment to either an entire regime (up-front) or a stage 1 treatment
(sequential) according to probabilities that are updated at $t$ based on the
accrued data from subjects previously enrolled at weeks $t-1,t-2,\ldots,0$.
Each of these subjects is assigned a regime (up-front) or $A_{1}$ (sequential)
using these probabilities. At each subsequent stage $k=2,\ldots,K$, we assume
that there is a delay, measured in weeks, between when $A_{k-1}$ is assigned
and when $\boldsymbol{X}_{k}$ is ascertained and $A_{k}$ is assigned, as well
as a follow-up period after $A_{K}$ is assigned but before $Y$ is recorded.
Thus, under sequential randomization, these subjects will require
randomization to $A_{k}$, $k=2,\ldots,K$, at a future week $(t+s)$, say, using
probabilities based on the accrued data from subjects previously enrolled at
or before week $(t+s-1)$.

To represent the data available on subjects who are already enrolled at any
week $t$, let $\Gamma_{t}\in\\{0,1\\}$ be the indicator that a subject has
enrolled in the SMART and been assigned to a regime or stage 1 treatment by
week $t$ (i.e., at a week $\leq t$). For such subjects, let $\tau\leq t$ be
the week of enrollment and assignment; $\kappa_{t}\in\\{1,\ldots,K\\}$ be the
most recent stage reached by week $t$; and $\Delta_{t}\in\\{0,1\\}$ be the
indicator that a subject has completed follow up by week $t$, so that $Y$ has
been observed. The data available on a subject at week $t$ are then
$D_{t}=\Gamma_{t}\\{1,\tau,\kappa_{t},\boldsymbol{X}_{1},A_{1},I(\kappa_{t}>1)\boldsymbol{X}_{2},I(\kappa_{t}>1)A_{2},\ldots,I(\kappa_{t}>K-1)\boldsymbol{X}_{K},I(\kappa_{t}>K-1)A_{K},\Delta_{t},\Delta_{t}Y\\}$;
$D_{t}$ is null if a subject has not yet enrolled by $t$. The accrued data
from all previously-enrolled subjects that can inform the randomization
probabilities at week $t$ are then $\mathfrak{D}_{t-1}=\\{D_{t-1,i},\forall
i\mbox{ such that }\Gamma_{t-1,i}=1\\}$.

###  2.3 Thompson sampling

The central idea of TS is to map a treatment’s randomization probability to
the “belief” that it is optimal among the available options. This belief is
represented conventionally by the posterior probability that a treatment is
optimal in the sense that it optimizes expected outcome (Thompson, 1933). For
example, in the up-front case (and similarly for a single-stage treatment),
let $\widehat{\rho}_{t}^{j}\in[0,1]$ be the estimated belief that treatment
regime $j=1,\ldots,m$ is optimal at week $t$ based on the accrued data
$\mathfrak{D}_{t-1}$ from previously-enrolled subjects. Let $r_{t}^{j}$ be
regime $j$’s randomization probability for subjects needing randomization at
week $t$. Ordinarily, $r_{t}^{j}$ is taken to be a monotone function of
$\widehat{\rho}_{t}^{j}$; a popular choice is
$r_{t}^{j}=(\widehat{\rho}_{t}^{j})^{c_{t}}/\sum_{v=1}^{m}(\widehat{\rho}_{t}^{v})^{c_{t}},$
where $c_{t}\in[0,1]$ is a damping constant (Thall and Wathen, 2007). Smaller
values of $c_{t}$ pull the probabilities toward uniform randomization; higher
values pull them closer to the beliefs. The $c_{t}$ can be the same for all
$t$ or increase with $t$ to impose greater adaptation as data accumulate.
Because aggressive adaptation can lead to randomization probabilities
approaching 0 or 1 for some regimes or treatments, limiting exploration of all
options, one may impose clipping constants, i.e., lower/upper bounds on the
probabilities (Zhang et al., 2021a). This practice is consistent with the
positivity assumption.

In a fully Bayesian formulation, beliefs $\widehat{\rho}_{t}^{j}$ are
estimated posterior probabilities based on $\mathfrak{D}_{t-1}$. We propose a
frequentist analog based on the so-called confidence distribution (Xie and
Singh, 2013). Here, we provide a basic overview of the approach; details for
the up-front and sequential algorithms are in subsequent sections. Let
$\boldsymbol{\theta}$ be the vector of parameters that, with a subject’s
current history, determines a subject’s optimal treatment; e.g., under up-
front randomization,
$\boldsymbol{\theta}=(\theta^{1},\ldots,\theta^{m})^{T}=\\{\mathcal{V}(\mathbf{d}^{1}),\ldots,\mathcal{V}(\mathbf{d}^{m})\\}^{T}$.
Let $\widehat{\mathcal{P}}_{\boldsymbol{\theta},t}$ be the estimated
confidence distribution for $\boldsymbol{\theta}$ at week $t$ based on
$\mathfrak{D}_{t-1}$, e.g., the estimated asymptotic distribution of some
estimator $\widehat{\boldsymbol{\theta}}_{t}$ for $\boldsymbol{\theta}$, and
let
$\widetilde{\boldsymbol{\theta}}^{b}_{t}=(\widetilde{\theta}^{1,b}_{t},\ldots,\widetilde{\theta}^{m,b}_{t})^{T}$,
$b=1,\ldots,B$, be independent draws from
$\widehat{\mathcal{P}}_{\boldsymbol{\theta},t}$. Then for $j=1,\ldots,m$,
$\widehat{\rho}_{t}^{j}=B^{-1}\sum_{b=1}^{B}I\\{\widetilde{\theta}^{j,b}_{t}=\max(\widetilde{\theta}^{1,b}_{t},\ldots,\widetilde{\theta}^{m,b}_{t})\\}$.

At the onset of a SMART, a burn-in period of nonadaptive randomization may be
required from which to obtain an initial estimate of $\boldsymbol{\theta}$ and
associated confidence distribution to be used for the first update of
randomization probabilities. The burn-in may be characterized in terms of
numbers of subjects who have completed each stage or the calendar time elapsed
since the start of the trial; examples are presented in Section 5. In what
follows, $t^{*}$ is the calendar time at which burn-in is complete, so that
adaptation starts at week $t=t^{*}+1$.

##  3 Adaptive randomization for SMARTs using Thompson sampling

We now present the proposed up-front and sequential TS approaches to RAR for
SMARTs. As above, under up-front randomization, subjects are randomized once,
at enrollment, to an embedded regime, which they follow through all $K$
stages. Thus, the “path” that a subject will take through the trial is
determined by $\mathfrak{D}_{t-1}$. Up-front randomization is logistically
simpler but does not use additional data that have accumulated as the subject
progresses through the trial. Up-front randomization is preferred when
simplicity of implementation is a priority or when enrollment is expected to
be slow relative to the time to progress through the trial, so that the amount
of new data accumulating before a subject completes all $K$ stages is modest.
Under sequential randomization, subjects are randomized at each stage, so that
up-to-date information from previous subjects informs randomization
probabilities as a subject progresses, but involves greater logistical
complexity.

###  3.1 Up-front randomization among regimes

To randomize newly-enrolled subjects at week $t$ to the embedded regimes
$\mathbf{d}^{j}$, $j=1,\ldots,m$, we require an estimator
$\widehat{\boldsymbol{\theta}}_{t}$ for
$\boldsymbol{\theta}=(\theta^{1},\ldots,\theta^{m})^{T}=\\{\mathcal{V}(\mathbf{d}^{1}),\ldots,\mathcal{V}(\mathbf{d}^{m})\\}^{T}$
based on $\mathfrak{D}_{t-1}$ with which to construct a confidence
distribution and thus randomization probabilities $r_{t}^{j}$. Aligned with
standard methods for SMARTs, we focus on inverse probability weighted (IPW)
and augmented IPW (AIPW) estimators for $\theta^{1},\ldots,\theta^{m}$
(Tsiatis et al., 2020, Sections 6.4.3-6.4.4). As we discuss in Section 4,
under any form of RAR, the asymptotic distribution of these estimators based
on the data at the end of the trial need not be normal; thus, we adapt the
approach of Zhang et al. (2021b) and propose weighted versions of these
estimators, where the weights are chosen so that the weighted estimators are
asymptotically normal. For stratified sampling, the formulation applies within
each stratum.

For subjects who enroll at week $t$, define
$\eta_{t,1}(a_{1},\boldsymbol{x}_{1},\mathfrak{D}_{t-1})=P(A_{1}=a_{1}|\boldsymbol{X}_{1}=\boldsymbol{x}_{1},\mathfrak{D}_{t-1})=P(A_{1}=a_{1}|\boldsymbol{H}_{1}=\boldsymbol{h}_{1},\mathfrak{D}_{t-1})$
and, for $k=2,\ldots,K$,
$\eta_{t,k}(a_{k},\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k-1},\mathfrak{D}_{t-1})=P(A_{k}=a_{k}|\overline{\boldsymbol{X}}_{k}=\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{A}}_{k-1}=\overline{\boldsymbol{a}}_{k-1},\mathfrak{D}_{t-1})=P(A_{k}=a_{k}|\boldsymbol{H}_{k}=\boldsymbol{h}_{k},\mathfrak{D}_{t-1})$.
For regime $\mathbf{d}$, let
$\overline{d}_{1}(\boldsymbol{x}_{1})=d_{1}(\boldsymbol{x}_{1})$, $\ldots$,
$\overline{d}_{k}(\overline{\boldsymbol{x}}_{k})=[d_{1}(\boldsymbol{x}_{1}),d_{2}\\{\overline{\boldsymbol{x}}_{2},d_{1}(\boldsymbol{x}_{1})\\}$,
$\ldots$,
$d_{k}\\{\overline{\boldsymbol{x}}_{k},\overline{d}_{k-1}(\overline{\boldsymbol{x}}_{k-1})\\}]$,
$k=3,\ldots,K$. For each $\mathbf{d}^{j}$, let
$C^{j}=I\\{\overline{\boldsymbol{A}}_{K}=\overline{d}_{K}^{j}(\overline{\boldsymbol{X}}_{K})\\}$
be the indicator that a subject’s experience through all $K$ stages is
consistent with receiving treatment using $\mathbf{d}^{j}$. Each enrolled
subject in $\mathfrak{D}_{t-1}$ was randomized at week $\tau\leq t-1$ based on
$\mathfrak{D}_{\tau-1}\subseteq\mathfrak{D}_{t-1}$, thus using $r_{\tau}^{j}$,
$j=1,\ldots,m$. For each, define
$\pi_{\tau,1}^{j}(\boldsymbol{x}_{1})=\eta_{\tau,1}\\{d_{1}^{j}(\boldsymbol{x}_{1}),\boldsymbol{x}_{1},\mathfrak{D}_{\tau-1}\\}$,
$\pi_{\tau,k}^{j}(\overline{\boldsymbol{x}}_{k})=\eta_{\tau,k}[d_{k}^{j}\\{\overline{\boldsymbol{x}}_{k},\overline{d}_{k-1}^{j}(\overline{\boldsymbol{x}}_{k-1})\\},\overline{\boldsymbol{x}}_{k},\overline{d}_{k-1}^{j}(\overline{\boldsymbol{x}}_{k-1}),\mathfrak{D}_{\tau-1}]$,
$k=2,\ldots,K$. The weighted IPW (WIPW) estimator for $\theta^{j}$ is

$\widehat{\theta}^{\mathrm{WIPW},j}_{t}=\left[\sum^{N}_{i=1}\frac{W_{\tau_{i}}^{j}\Delta_{t-1,i}C_{i}^{j}}{\\{\prod_{k=2}^{K}\pi_{\tau_{i},k}^{j}(\overline{\boldsymbol{X}}_{k,i})\\}\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})}\right]^{-1}\\!\\!\\!\sum^{N}_{i=1}\frac{W_{\tau_{i}}^{j}\Delta_{t-1,i}C_{i}^{j}Y_{i}}{\\{\prod_{k=2}^{K}\pi_{\tau_{i},k}^{j}(\overline{\boldsymbol{X}}_{k,i})\\}\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})},$
(1)

where $W_{\tau}^{j}$ is a weight depending on $\mathfrak{D}_{\tau-1}$
discussed further in Section 4. The denominator in each term in (1) can be
interpreted as the propensity for receiving treatment consistent with
$\mathbf{d}^{j}$ through all $K$ stages and depends on $r_{\tau}^{j}$,
$j=1,\ldots,m$. For example, in the cancer pain SMART with $K=2$, let $X_{1}$
denote a patient’s baseline pain score and $X_{2,1}$ their pain score at the
end of stage 1, so that $\overline{\boldsymbol{X}}_{2}=(X_{1},X_{2,1})^{T}$,
and let $X_{2,2}=I(X_{2,1}\leq 0.7\,X_{1})$ denote response status. For
embedded regime $\mathbf{d}^{1}$, under which a subject will receive $A_{1}=0$
and then $A_{2}=0$ if they respond to $A_{1}$ and $A_{2}=1$ otherwise,
$C^{1}=I(A_{1}=0,X_{2,2}=1,A_{2}=0)+I(A_{1}=0,X_{2,2}=0,A_{2}=1)$. Because
regimes 1–4 assign $A_{1}=0$,
$\pi_{\tau,1}^{1}(\boldsymbol{X}_{1})=\sum^{4}_{j=1}r_{\tau}^{j}$, and because
regimes 1 and 2 assign $A_{2}=0$ to responders and regimes 1 and 4 assign
$A_{2}=1$ to nonresponders,
$\pi_{\tau,2}^{1}(\overline{\boldsymbol{X}}_{2})=\\{I(X_{2,2}=1)(r_{\tau}^{1}+r_{\tau}^{2})+I(X_{2,2}=0)(r_{\tau}^{1}+r_{\tau}^{4})\\}/\sum^{4}_{j=1}r_{\tau}^{j}$.

Usual AIPW estimators incorporate baseline and interim information (e.g.,
Tsiatis et al., 2020, Section 6.4.4) to gain efficiency over IPW estimators.
Accordingly, we consider a class of weighted AIPW (WAIPW) estimators. Let
$W_{\tau}^{A,j}$ be a weight depending on $\mathfrak{D}_{\tau-1}$; and
$\overline{C}_{k}^{j}=I\\{\overline{\boldsymbol{A}}_{k}=\overline{d}_{k}^{j}(\overline{\boldsymbol{X}}_{k})\\}$,
$k=1,\ldots,K$, with $\overline{C}_{0}^{j}\equiv 1$. Estimators in the class
are of the form

$\displaystyle\widehat{\theta}_{t}^{\mathrm{WAIPW},j}=\left(\sum^{N}_{i=1}W_{\tau_{i}}^{A,j}\Delta_{t-1,i}\right)^{-1}\sum^{N}_{i=1}W_{\tau_{i}}^{A,j}\Delta_{t-1,i}\left(\frac{C_{i}^{j}Y_{i}}{\left\\{\prod_{k=2}^{K}\pi_{\tau_{i},k}^{j}(\overline{\boldsymbol{X}}_{k,i})\right\\}\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})}\vphantom{\frac{\overline{C}_{k,i}^{j}}{\left\\{\prod_{v=2}^{k-1}\pi_{d,v,tau_{i}}(\overline{\boldsymbol{X}}_{v,i})\right\\}\pi_{d,1\tau_{i}}(\boldsymbol{X}_{1,i})}}\right.$
(2)
$\displaystyle\,\,+\left.\sum_{k=1}^{K}\left[\frac{\overline{C}_{k-1,i}^{j}}{\left\\{\prod_{v=2}^{k-1}\pi_{\tau_{i},v}^{j}(\overline{\boldsymbol{X}}_{v,i})\right\\}\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})}-\frac{\overline{C}_{k,i}^{j}}{\left\\{\prod_{v=2}^{k}\pi_{\tau_{i},v}^{j}(\overline{\boldsymbol{X}}_{v,i})\right\\}\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})}\right]L^{j}_{k}(\overline{\boldsymbol{X}}_{k,i})\right),$

where $L^{j}_{k}(\overline{\boldsymbol{x}}_{k})$ is an arbitrary function of
$\overline{\boldsymbol{x}}_{k}$ and $\mathbf{d}^{j}$, $k=1,\ldots,K$; and
$\prod_{v=1}^{0}\pi_{\tau,v}^{j}(\boldsymbol{x}_{v})\equiv 1$.

The optimal choice is
$L^{j}_{k}(\overline{\boldsymbol{x}}_{k})=E\\{Y^{*}(\mathbf{d}^{j})|\overline{\boldsymbol{X}}_{k}=\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{A}}_{k-1}=\overline{d}_{k-1}^{j}(\overline{\boldsymbol{x}}_{k-1})\\}$,
$k=1,\ldots,K$, which can be modeled and estimated by adapting the backward
iterative scheme in Tsiatis et al. (2020, Section 6.4.2). Just as the
denominators of each term in (1) and (2) depend on $r_{\tau}^{j}$,
$j=1,\ldots,m$, based on $\mathfrak{D}_{\tau-1}\subseteq\mathfrak{D}_{t-1}$,
the fitted models used to approximate the optimal
$L^{j}_{k}(\overline{\boldsymbol{X}}_{k})$ should be based on
$\mathfrak{D}_{\tau-1}$. At stage $K$, define
$Q_{K}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K})=E(Y|\overline{\boldsymbol{X}}_{K}=\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{A}}_{K}=\overline{\boldsymbol{a}}_{K})$
and
$V_{K}^{j}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K-1})=Q_{K}\\{\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K-1},d_{K}^{j}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K-1})\\}$.
Posit a model
$Q_{K}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K};\boldsymbol{\beta}_{K})$
for $Q_{K}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K})$
indexed by $\boldsymbol{\beta}_{K}$, e.g., a linear or logistic model for
continuous or binary $Y$, respectively. Let
$\widehat{\boldsymbol{\beta}}_{K,t}$ denote an estimator for
$\boldsymbol{\beta}_{K}$ based on subjects in $\mathfrak{D}_{t-1}$ for whom
$\Delta_{t-1}=1$, e.g., using least squares or maximum likelihood, and define
$\widehat{Q}_{K,t}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K})=Q_{K}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K};\widehat{\boldsymbol{\beta}}_{K,t})$
and
$\widetilde{V}_{K,t}^{j}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K-1})=\widehat{Q}_{K,t}\\{\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K-1},d_{K}^{j}(a_{K})\\}$.
Recursively for $k=K-1,\ldots,1$, define
$Q_{k}^{j}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k})=\allowbreak
E\\{V_{k+1}^{j}(\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k})|\overline{\boldsymbol{X}}_{k}=\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{A}}_{k}=\overline{\boldsymbol{a}}_{k}\\}$
and
$V_{k}^{j}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k-1})=Q_{k}^{j}\\{\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k-1},d_{k}^{j}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k-1})\\}$.
Posit a model
$Q_{k}^{j}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k};\boldsymbol{\beta}_{k}^{j})$
indexed by $\boldsymbol{\beta}_{k}^{j}$, and obtain estimator
$\widehat{\boldsymbol{\beta}}_{k,t}^{j}$ based on subjects in
$\mathfrak{D}_{t-1}$ for whom $\Delta_{t-1}=1$ by an appropriate regression
method using as the outcome the pseudo outcomes
$\widetilde{V}_{k+1,t}^{j}(\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k})=\widehat{Q}_{k+1,t}^{j}\\{\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k},d_{k+1}^{j}(\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k})\\}$,
and let
$\widehat{Q}_{k,t}^{j}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k})=Q_{k}^{j}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k};\widehat{\boldsymbol{\beta}}_{k,t}^{j})$.
As in Tsiatis et al. (2020, Section 6.4.2), these models may involve separate
expressions for responders and nonresponders; and, if at stage $k$ there is
only one treatment option for a subject’s history, $Y$ (if $k=K-1$) or
$\widetilde{V}^{j}_{k+2,t}(\overline{\boldsymbol{X}}_{k+2},\overline{\boldsymbol{A}}_{k+1})$
(if $k<K-1$) can be “carried back” in place of
$\widetilde{V}_{k+1,t}^{j}(\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k})$.
For each $i$ with $\Delta_{t-1,i}=1$ in (2), substitute
$Q_{K}\\{\overline{\boldsymbol{X}}_{K,i},\overline{d}_{K}^{j}(\overline{\boldsymbol{X}}_{K,i});\widehat{\boldsymbol{\beta}}_{K,\tau_{i}}\\}$
and
$Q^{j}_{k}\\{\overline{\boldsymbol{X}}_{k,i},\overline{d}_{k}^{j}(\overline{\boldsymbol{X}}_{k,i});\widehat{\boldsymbol{\beta}}^{j}_{k,\tau_{i}}\\}$
for $L^{j}_{K}(\overline{\boldsymbol{X}}_{K,i})$ and
$L^{j}_{k}(\overline{\boldsymbol{X}}_{k,i})$, $k=1,\ldots,K-1$, respectively.

As the basis for RAR, we propose using (1) or (2), with or without weights, to
obtain estimators $\widehat{\theta}_{t}^{j}$ for
$\theta^{j}=\mathcal{V}(\mathbf{d}^{j})$, $j=1,\ldots,m$, based on
$\mathfrak{D}_{t-1}$; and take the estimated confidence distribution for
$\boldsymbol{\theta}$ needed to obtain $\widehat{\rho}_{t}^{j}$ to form
$r_{t}^{j}$, $j=1,\ldots,m$, to be the asymptotic normal distribution for
$\widehat{\boldsymbol{\theta}}_{t}=\\{\widehat{\theta}_{t}^{1},\ldots,\widehat{\theta}_{t}^{m}\\}^{T}$
following from M-estimation theory (e.g., Tsiatis et al., 2020, Section
6.4.4), with the weights treated as fixed.

Basing confidence distributions and thus randomization probabilities at week
$t$ on (1) or (2) uses data only on subjects in $\mathfrak{D}_{t-1}$ who have
completed the trial, $\Delta_{t-1}=1$. To exploit partial information on
subjects still progressing through the trial at $t$, with $\Delta_{t-1}=0$, it
is possible to develop a weighted version of the interim AIPW (IAIPW)
estimator of Manschot, Laber, and Davidian (2023); see Appendix A. Simulations
in Section 5 show negligible gains in performance over (1) or (2).

###  3.2 Sequential randomization based on the optimal regime

We propose methods for obtaining randomization probabilities at week $t$ based
on $\mathfrak{D}_{t-1}$ to be used to assign treatments for subjects requiring
randomization at $t$ at any stage $k=1,\ldots,K$. Because the set of feasible
treatments $\Psi_{k}(\boldsymbol{h}_{k})$ for a subject with history
$\boldsymbol{h}_{k}$ at stage $k$ may depend on $\boldsymbol{h}_{k}$, as when
the sets of options for responders and nonresponders to previous treatment are
different, randomization probabilities may be history dependent. The approach
uses Q-learning for estimation of an optimal, individualized regime (e.g.,
Tsiatis et al., 2020, Section 7.4.1). We present the approach when $Y$ is
continuous and linear models are used; extensions to other outcomes and more
flexible models are possible Moodie et al. (2014).

Define the Q-functions
$Q_{K}(\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{a}}_{K})=E(Y\mid\overline{\boldsymbol{X}}_{K}=\overline{\boldsymbol{x}}_{K},\overline{\boldsymbol{A}}_{K}=\overline{\boldsymbol{a}}_{K})$
and, for $k=K-1,\ldots,1$,
$Q_{k}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k})=E\left\\{V_{k+1}(\overline{\boldsymbol{x}}_{k},\boldsymbol{X}_{k+1},\overline{\boldsymbol{a}}_{k})\mid\overline{\boldsymbol{X}}_{k}=\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{A}}_{k}=\overline{\boldsymbol{a}}_{k}\right\\}$,
where
$V_{k}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k-1})=\max_{a_{k}\in\Psi_{k}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k-1})}Q_{k}(\overline{\boldsymbol{x}}_{k},\overline{a}_{k-1},a_{k})$,
$k=1,\ldots,K$. Posit models
$Q_{k}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k};\boldsymbol{\beta}_{k})=\boldsymbol{\phi}_{k}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k})^{T}\boldsymbol{\beta}_{k}$,
where
$\boldsymbol{\phi}_{k}(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k})$
is a $p_{k}$-dimensional feature function, $k=1,\ldots,K$. Randomization
probabilities are obtained via the following backward algorithm. At stage $K$,
obtain $\widehat{\boldsymbol{\beta}}_{K,t}$, the ordinary least squares (OLS)
estimator based on subjects in $\mathfrak{D}_{t-1}$ for whom $\Delta_{t-1}=1$
and its estimated covariance matrix
$\widehat{\boldsymbol{\Sigma}}_{K,t}=\widehat{\sigma}^{2}_{t,K}\big{\\{}\sum^{N}_{i=1}\Delta_{t-1,i}\boldsymbol{\phi}_{K}(\overline{\boldsymbol{X}}_{K,i},\overline{\boldsymbol{A}}_{K,i})\boldsymbol{\phi}_{K}(\overline{\boldsymbol{X}}_{K,i},\overline{\boldsymbol{A}}_{K,i})^{T}\big{\\}}^{-1}$,
$\widehat{\sigma}^{2}_{K,t}=N_{t}^{-1}\sum^{N}_{i=1}\Delta_{t-1,i}\\{Y_{i}-\boldsymbol{\phi}_{K}(\overline{\boldsymbol{X}}_{K,i},\overline{\boldsymbol{A}}_{K,i})^{T}\widehat{\boldsymbol{\beta}}_{K,t}\\}^{2}$,
$N_{t}=\sum^{N}_{i=1}\Delta_{t-1,i}$. Based on these results, obtain the
estimated confidence distribution
$\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{K,t}}$ for
$\boldsymbol{\beta}_{K}$ as described below. For $k=K-1,\ldots,1$, define
pseudo outcomes
$\widetilde{V}_{k+1,t}(\boldsymbol{\beta}_{k+1})=\max_{a_{k+1}\in\Psi(\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k})}Q_{k+1}(\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k},a_{k+1};\boldsymbol{\beta}_{k+1})$
for subjects in $\mathfrak{D}_{t-1}$ with $\Delta_{t-1}=1$. If a subject’s
history is such that
$\Psi_{k+1}(\overline{\boldsymbol{X}}_{k+1},\overline{\boldsymbol{A}}_{k})$
comprises a single treatment option,
$\widetilde{V}_{k+1,t}(\boldsymbol{\beta}_{k+1})$ can be taken equal to the
pseudo outcome at step $k+2$ or $Y$ if $k=K-1$. Obtain an estimator for
$\boldsymbol{\beta}_{k}$, $k=K-1,\ldots,1$, by OLS as

$\widehat{\boldsymbol{\beta}}_{k,t}(\widehat{\boldsymbol{\beta}}_{k+1,t})=\mbox{argmin}_{\boldsymbol{\beta}}\sum^{N}_{i=1}\Delta_{t-1,i}\\{\widetilde{V}_{k+1,t,i}(\widehat{\boldsymbol{\beta}}_{k+1,t})-\boldsymbol{\phi}_{k}(\overline{\boldsymbol{X}}_{k,i},\overline{\boldsymbol{A}}_{k,i})^{T}\boldsymbol{\beta}\\}^{2},$
(3)

and estimated covariance matrix
$\widehat{\boldsymbol{\Sigma}}_{k,t}(\widehat{\boldsymbol{\beta}}_{k+1,t})=\widehat{\sigma}^{2}_{k,t}\big{\\{}\sum^{N}_{i=1}\Delta_{t-1,i}\boldsymbol{\phi}_{k}(\overline{\boldsymbol{X}}_{k,i},\overline{\boldsymbol{A}}_{k,i})\boldsymbol{\phi}_{k}(\overline{\boldsymbol{X}}_{k,i},\overline{\boldsymbol{A}}_{k,i})^{T}\big{\\}}^{-1}$,
$\widehat{\sigma}^{2}_{k,t}(\widehat{\boldsymbol{\beta}}_{k+1,t})=N_{t}^{-1}\sum^{N}_{i=1}\Delta_{t-1,i}\\{\widetilde{V}_{k+1,t,i}(\widehat{\boldsymbol{\beta}}_{k+1,t})-\boldsymbol{\phi}_{k}(\overline{\boldsymbol{X}}_{k,i},\overline{\boldsymbol{A}}_{k,i})^{T}\widehat{\boldsymbol{\beta}}_{k,t}(\widehat{\boldsymbol{\beta}}_{k+1,t})\\}^{2}$.
Based on these results, obtain an estimated confidence distribution for
$\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{k,t}}$ as described next. Note
that, if subjects enter the SMART by a random process, basing the fitted
models on subjects in $\mathfrak{D}_{t-1}$ with $\Delta_{t-1}$ is reasonable,
as these subjects are representative of the subject population.

Because $\widehat{\boldsymbol{\beta}}_{K,t}$ is a standard OLS estimator, it
is natural to approximate the confidence distribution
$\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{K,t}}$ for
$\boldsymbol{\beta}_{K}$ by
$\mathcal{N}(\widehat{\boldsymbol{\beta}}_{K,t},\widehat{\boldsymbol{\Sigma}}_{K,t})$.
However, because (3) is not a standard regression problem,
$\widehat{\boldsymbol{\beta}}_{k,t}(\widehat{\boldsymbol{\beta}}_{k+1,t})$,
$k<K$, is a nonregular estimator (Tsiatis et al., 2020, Section 10.4.1). Thus,
usual large sample theory does not apply and confidence intervals based on
(unadjusted) normal or bootstrap approximations need not achieve nominal
coverage. Accordingly, we obtain a confidence distribution
$\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{k,t}}$ for
$\boldsymbol{\beta}_{k}$, $k=1,\ldots,K-1$, in the spirit of a projection
interval (Laber et al., 2014), which faithfully represents the uncertainty in
$\widehat{\boldsymbol{\beta}}_{K-1,t},\dots,\widehat{\boldsymbol{\beta}}_{1,t}$.

We demonstrate this approach for $K=3$. For final stage 3, draw a sample of
size $B_{3}$,
$\widetilde{\boldsymbol{\beta}}_{3,t}^{1},\ldots,\widetilde{\boldsymbol{\beta}}_{3,t}^{B_{3}}$,
say, from $\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{3,t}}$, the approximate
normal sampling distribution
$\mathcal{N}(\widehat{\boldsymbol{\beta}}_{3,t},\widehat{\boldsymbol{\Sigma}}_{3,t})$
as above. At stage 2, first draw a sample
$\widetilde{\boldsymbol{\beta}}_{3,t}^{1},\ldots,\widetilde{\boldsymbol{\beta}}_{3,t}^{b_{3}}$
of size $b_{3}$ from $\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{3},t}$. For
each $\widetilde{\boldsymbol{\beta}}_{3,t}^{b}$, form pseudo outcomes
$\widetilde{V}_{3,t,i}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b})=\max_{a_{3}\in\Psi(\overline{\boldsymbol{X}}_{3,i},\overline{\boldsymbol{A}}_{2,i})}Q_{3}(\overline{\boldsymbol{X}}_{3,i},\overline{\boldsymbol{A}}_{2,i},a_{3};\widetilde{\boldsymbol{\beta}}_{3,t}^{b})$
and obtain
$\widehat{\boldsymbol{\beta}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b})$
by OLS analogous to (3) with $k=2$, and obtain
$\widehat{\boldsymbol{\Sigma}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b})$
similarly. Then draw a sample
$\widetilde{\boldsymbol{\beta}}_{2,t}^{1},\ldots,\widetilde{\boldsymbol{\beta}}_{2,t}^{b_{2}}$
from
$\mathcal{N}\\{\widehat{\boldsymbol{\beta}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b}),\widehat{\boldsymbol{\Sigma}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b})\\}$.
The $b_{3}$ samples of size $b_{2}$ corresponding to each
$\widetilde{\boldsymbol{\beta}}_{3,t}^{b}$, $b=1,\ldots,b_{3}$, collectively
comprise a sample of size $B_{2}=b_{3}\times b_{2}$ from the confidence
distribution $\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{2,t}}$ for
$\boldsymbol{\beta}_{2}$. At stage 1, again draw a sample
$\widetilde{\boldsymbol{\beta}}_{3,t}^{1},\ldots,\widetilde{\boldsymbol{\beta}}_{3,t}^{b_{3}}$
from $\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{3},t}$ and obtain
$\widehat{\boldsymbol{\beta}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b})$
and
$\widehat{\boldsymbol{\Sigma}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b})$,
$b=1,\ldots,b_{3}$, as above. Then draw a sample
$\widetilde{\boldsymbol{\beta}}_{2,t}^{1},\ldots,\widetilde{\boldsymbol{\beta}}_{2,t}^{b_{2}}$
from
$\mathcal{N}\\{\widehat{\boldsymbol{\beta}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b}),\widehat{\boldsymbol{\Sigma}}_{2,t}(\widetilde{\boldsymbol{\beta}}_{3,t}^{b})\\}$.
For each of $\widetilde{\boldsymbol{\beta}}_{2,t}^{b}$, $b=1,\ldots,b_{2}$,
form pseudo outcomes
$\widetilde{V}_{2,t,i}(\widetilde{\boldsymbol{\beta}}_{2,t}^{b})=\max_{a_{2}\in\Psi(\overline{\boldsymbol{X}}_{2,i},A_{1,i})}Q_{2}(\overline{\boldsymbol{X}}_{2,i},A_{1,i},a_{2};\widetilde{\boldsymbol{\beta}}_{2,t}^{b})$,
and obtain
$\widehat{\boldsymbol{\beta}}_{1,t}(\widetilde{\boldsymbol{\beta}}_{2,t}^{b})$
by OLS analogous to (3) with $k=1$, and obtain
$\widehat{\boldsymbol{\Sigma}}_{t,1}(\widetilde{\boldsymbol{\beta}}_{2,t}^{b})$.
Then draw a sample
$\widetilde{\boldsymbol{\beta}}_{1,t}^{1},\ldots,\widetilde{\boldsymbol{\beta}}_{1,t}^{b_{1}}$
from
$\mathcal{N}\\{\widehat{\boldsymbol{\beta}}_{1,t}(\widetilde{\boldsymbol{\beta}}_{2,t}^{b}),\widehat{\boldsymbol{\Sigma}}_{1,t}(\widetilde{\boldsymbol{\beta}}_{2,t}^{b})\\}$,
$b=1,\ldots,b_{2}$. The $b_{3}\times b_{2}$ samples of size $b_{1}$
corresponding to each combination of draws
$\widetilde{\boldsymbol{\beta}}_{3,t}^{b}$,
$\widetilde{\boldsymbol{\beta}}_{2,t}^{b^{\prime}}$, $b=1,\ldots,b_{3}$,
$b^{\prime}=1,\ldots,b_{2}$, collectively comprise a sample of size
$B_{1}=b_{3}\times b_{2}\times b_{1}$ from the confidence distribution
$\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{1},t}$. This procedure is
embarrassingly parallel, and as it involves only draws from a multivariate
normal distribution, it is not computationally burdensome.

Having generated $B_{k}$ draws from
$\widehat{\mathcal{P}}_{\boldsymbol{\beta}_{t,k}}$, $k=1,\ldots,K$, at week
$t$, we can obtain beliefs, i.e., draws from approximate posterior
distributions, which can be translated into randomization probabilities for
subjects requiring treatment assignments at any stage $k=1,\ldots,K$ at week
$t$. Because the $Q$-functions depend on patient history, the randomization
probabilities under TS can vary across subjects at a given time even if their
feasible treatments are the same. To see this, suppose that there are $s_{k}$
feasible sets of treatments at stage $k$, denoted by
$\zeta_{k}^{u}=\\{a_{k}^{1,u},\ldots,a_{k}^{m_{k}^{u},u}\\}$,
$u=1,\ldots,s_{k}$. Further, suppose that there are $\ell_{k,t}$ subjects
requiring randomization at stage $k$ at week $t$ and that $\ell_{k,t}^{u}$ are
eligible for the $u$th feasible set. For the $u$th feasible set, randomization
probabilities can be obtained for each subject $v=1,\ldots,\ell_{k,t}^{u}$ by
defining the belief for $a_{k}^{j,u}$, $j=1,\ldots,m_{k}^{u}$, depending on
$\mathfrak{D}_{t-1}$ for subject $v$ as
$\widehat{\rho}_{t,v}^{j,k,u}=B_{k}^{-1}\sum^{B_{k}}_{b=1}I\\{a_{k}^{j,u}=\mbox{argmax}_{a\in\zeta_{k}^{u}}Q_{k}(\overline{\boldsymbol{X}}_{k,v},\overline{\boldsymbol{A}}_{k-1,v},a;\widetilde{\boldsymbol{\beta}}_{k,t}^{b})\\},$
from which randomization probabilities for subject $v$ are obtained. If the
Q-function model depends on the history
$\boldsymbol{h}_{k}=(\overline{\boldsymbol{x}}_{k},\overline{\boldsymbol{a}}_{k-1})$
only through the components such as previous treatment and response status
that dictate the feasible set, then this approach will yield the same
probabilities for all $\ell_{k,t}^{u}$ subjects. Otherwise, randomization
probabilities will be individual-specific, depending on covariate and
treatment information in addition to $\mathfrak{D}_{t-1}$, which could be
logistically complex. A second approach is to consider all
$(m_{k}^{u})^{\ell_{k,t}^{u}}$ possible configurations for assigning the
options in $\zeta_{k}^{u}$ to the $\ell_{k,t}^{u}$ subjects and define the
beliefs and thus randomization probabilities based on the configuration that
maximizes the average $k$th Q-function across subjects. Letting
$\Lambda_{k}^{u}$ be the set of all possible configurations, writing the $j$th
configuration as $\mbox{\sf a}^{j}=(\mbox{\sf a}_{1}^{j},\ldots,\mbox{\sf
a}_{\ell_{k,t}^{u}}^{j})$, define
$\widehat{\rho}_{t}^{j,k,u}=B_{k}^{-1}\sum_{b=1}^{B_{k}}I\Big{\\{}\mbox{\sf
a}^{j}=\mbox{argmax}_{\mbox{\sf
a}\in\Lambda_{k}^{u}}(\ell_{k,t}^{u})^{-1}\sum_{v=1}^{\ell_{k,t}^{u}}Q_{k}(\overline{\boldsymbol{X}}_{k,v},\overline{\boldsymbol{A}}_{k-1,v},\mbox{\sf
a};\widetilde{\boldsymbol{\beta}}_{k,t}^{b})\Big{\\}}.$ Simulation experiments
(not shown here) suggest that the two approaches perform similarly.

##  4 Post-Trial Inference

Although the potential outcomes $\mathcal{W}^{*}_{i}$, $i=1,\ldots,N$, are
independent and identically distributed (i.i.d.), under any form of RAR, the
observed data are not, as the randomization probabilities are functions of the
past data $\mathfrak{D}_{t-1}$ (Zhang et al., 2021b). Thus, post-trial
evaluation of $\mathcal{V}(\mathbf{d}^{1}),\ldots,\mathcal{V}(\mathbf{d}^{m})$
based on the usual unweighted IPW or AIPW estimators is potentially
problematic, as standard asymptotic theory for these estimators, which assumes
i.i.d. data, does not apply (e.g., Bibaut et al., 2021; Zhang et al., 2021a).
Thus, we adapt the approach of Zhang et al. (2021b) and choose the weights in
(1) and (2) so that asymptotic normality for these estimators can be
established via the martingale central limit theorem We sketch the rationale;
details are given in Appendix B.

To emphasize the key ideas, consider a simplified setting in which a fixed
number of subjects, $n$, enrolls at each week and $T\rightarrow\infty$, so
that the total number of subjects $N=nT\rightarrow\infty$. To simplify
notation, take $n=1$. At the end of the trial, with all data complete,
reindexing subjects by $t=1,\ldots,T$, for regime $j$, the estimators
$\widehat{\theta}_{T}^{\mathrm{WIPW},j}$ and
$\widehat{\theta}_{T}^{\mathrm{WAIPW},j}$, say, are solutions in $\theta^{j}$
to an estimating equation
$\sum_{t=1}^{T}M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})=0$;
e.g., for (1),

$M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})=\frac{W_{t}^{j}C^{j}_{t}}{\\{\prod_{k=2}^{K}\pi_{t,k}^{j}(\overline{\boldsymbol{X}}_{k,t})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1,t})}(Y_{t}-\theta^{j}).$
(4)

Critical to the proof is that the weights $W_{t}^{j}$ and $W^{A,j}_{t}$ are
chosen so that (i) the estimating equations remain conditionally unbiased,
$E\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})|\mathfrak{D}_{t-1}\\}=0$,
and (ii) the variance is stabilized,
$E[\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})\\}^{2}|\mathfrak{D}_{t-1}]$
$=\sigma^{2}>0$ for all $t$. Writing $\widehat{\theta}_{T}^{j}$ to denote
either estimator and defining
$\delta_{T}^{j}=T^{-1}\sum_{t=1}^{T}\partial/\partial\theta^{j}\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y;\theta^{j})\\}_{\theta^{j}=\widehat{\theta}_{T}^{j}}$
and
$(\sigma^{j}_{T})^{2}=T^{-1}\sum_{t=1}^{T}M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\widehat{\theta}_{T}^{j})^{2}$,
we argue in Appendix B that, under standard regularity conditions,
$\widehat{\theta}_{T,j}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\theta^{j}$ and

$\delta^{j}_{T}(\sigma_{T}^{j})^{-1}T^{1/2}(\widehat{\theta}_{T}^{j}-\theta^{j})\stackrel{{\scriptstyle
D}}{{\longrightarrow}}\mathcal{N}(0,1).$ (5)

In Section 5, we use (5) to construct confidence intervals and bounds for
$\theta^{j}$, $j=1,\ldots,m$.

We sketch arguments for the WIPW estimator to show that (i) holds and how to
choose the weights to guarantee (ii); see Appendix B for details and arguments
for the WAIPW estimator. Define
$\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})=\\{\boldsymbol{X}_{1},\boldsymbol{X}_{2}^{*}(\mathbf{d}^{j}),\ldots,\boldsymbol{X}_{k}^{*}(\mathbf{d}^{j})\\}$,
$k=1,\ldots,K$, and write
$\pi_{t}^{*,j}\\{\overline{\boldsymbol{X}}_{K}^{*}(\mathbf{d}^{j})\\}=\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})$.
When $C^{j}_{t}=1$, $Y_{t}=Y^{*}(\mathbf{d}^{j})$,
$\overline{\boldsymbol{X}}_{K,t}=\overline{\boldsymbol{X}}_{K}^{*}(\mathbf{d}^{j})$,
so that (i) is

$\displaystyle E$
$\displaystyle\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})|\mathfrak{D}_{t-1}\\}=W_{t}^{j}E\left[\left.\frac{C^{j}_{t}}{\pi_{t}^{*,j}\\{\overline{\boldsymbol{X}}_{K}^{*}(\mathbf{d}^{j})\\}}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}\right|\mathfrak{D}_{t-1}\right]$
$\displaystyle=W_{t}^{j}E\left[\left.\frac{E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})}{\pi_{t}^{*,j}\\{\overline{\boldsymbol{X}}_{K}^{*}(\mathbf{d}^{j})\\}}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}\right|\mathfrak{D}_{t-1}\right]=W_{t}^{j}E\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}=0,$

because, as in Tsiatis et al. (2020, Section 6.4.3),
$E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\pi_{\tau}^{*,j}\\{\overline{\boldsymbol{X}}_{K}^{*}(\mathbf{d}^{j})\\}$,
and $\mathcal{W}^{*}$ and thus $Y^{*}(\mathbf{d}^{j})$ is independent of
$\mathfrak{D}_{t-1}$. Using similar manipulations,

$E\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})^{2}|\mathfrak{D}_{t-1}\\}=(W_{t}^{j})^{2}E\left[\left.\frac{\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}}{\pi_{t}^{*,j}\\{\overline{\boldsymbol{X}}_{K}^{*}(\mathbf{d}^{j})\\}}\right|\mathfrak{D}_{t-1}\right].$
(6)

Thus, to ensure (ii), $W_{t}^{j}$ should be chosen so that (6) is a constant
depending only on $j$.

We demonstrate the choice of $W_{t}^{j}$ in practice for regime 1 of the
cancer pain SMART, $K=2$, and up-front randomization. From Section 3.1,
letting $X^{*}_{2,2}(\mathbf{d}^{1})=I\\{X^{*}_{2,1}(\mathbf{d}^{1})\leq
0.7\,X_{1}\\}$,
$\pi_{t,1}^{1}(\boldsymbol{X}_{1})=\pi_{t,1}^{1}=\sum^{4}_{j=1}r_{t}^{j}$, and
$\pi_{t,2}^{1}\\{\overline{\boldsymbol{X}}^{*}_{2}(\mathbf{d}^{1})\\}=I\\{X_{2,2}^{*}(\mathbf{d}^{1})=1\\}\pi^{1(1)}_{t,2}+I\\{X_{2,2}^{*}(\mathbf{d}^{1})=0\\}\pi^{1(0)}_{t,2}$,
where $\pi^{1(1)}_{t,2}=(r_{t}^{1}+r_{t}^{2})/\pi_{t,1}^{1}$ and
$\pi^{1(0)}_{t,2}(r_{t}^{1}+r_{t}^{4})/\pi_{t,1}^{1}$. Then, using
$\mathcal{W}^{*}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathfrak{D}_{t-1}$, (6) is

$(W_{t,1})^{2}\left(\frac{\mu^{1(1)}}{\pi^{1(1)}_{t,2}\pi_{t,1}^{1}}+\frac{\mu^{(0)}_{1}}{\pi^{1(0)}_{t,2}\pi_{t,1}^{1}}\right),\,\,\,\mu^{1(s)}=E[I\\{X_{2,2}^{*}(\mathbf{d}^{1})=s\\}\\{Y^{*}(\mathbf{d}^{1})-\theta^{1}\\}^{2}],\,s=0,1.$
(7)

Then, in the original notation, estimate $\mu^{1(s)}_{1}$ at week $t$ by
$\widehat{\mu}^{1(s)}_{t}=N^{-1}_{t}\sum^{N}_{i=1}\big{[}\Delta_{t-1,i}I(X_{2,2,i}=s)C^{1}_{i}(Y_{i}-\widetilde{\theta}_{t}^{1})^{2}/\\{\pi^{1(s)}_{\tau_{i},2}\pi_{\tau_{i},1}^{1}\\}\big{]}$,
$s=0,1$, based on $\mathfrak{D}_{t-1}$, where $\widetilde{\theta}_{t}^{1}$ is
an estimator for $\theta^{1}$ using $\mathfrak{D}_{t-1}$ (we use the
unweighted IPW estimator). Setting (7) equal to a constant $\Xi^{1}$, say, and
defining
$\widehat{\Xi}_{t}^{1}=\sum_{s=0}^{1}\\{\widehat{\mu}^{1(s)}_{t}/(\pi^{1(s)}_{t,2}\pi_{t,1}^{1})\\}$,
estimate $\Xi^{1}$ based on the burn-in data, for which $W_{t}^{j}\equiv 1$
and $r_{t}^{j}$, $j=1,\ldots,8$, are fixed (nonadaptive) for $t\leq t^{*}$, by
$\widehat{\Xi}_{t^{*}}^{1}$. Then for $t\geq t^{*}+1$, take
$W_{t}^{1}=(\widehat{\Xi}^{1}_{t^{*}}/\widehat{\Xi}^{1}_{t})^{1/2}$. See
Appendix B for considerations for sequential RAR.

##  5 Simulation studies

###  5.1 Up-front randomization among embedded regimes

We present results of a simulation study involving 5000 Monte Carlo trials
under a scenario mimicking the cancer pain SMART introduced in Section 1 in
which subjects are randomized using various forms of up-front RAR as in
Section 3.1. Each trial enrolls $N=1000$ subjects, with enrollment times
uniform over (integer) weeks 1-24. Upon enrollment, we draw baseline pain
score $X_{1}\sim\mathcal{N}(5,1)$ and assign stage 1 treatment
$A_{1}\in\mathcal{A}_{1}=\\{0,1\\}$. Six weeks after $A_{1}$ is assigned,
second-stage pain score is generated as
$X_{2,1}=\gamma_{1,0}+\gamma_{1,1}X_{1}+\gamma_{1,2}A_{1}+\varepsilon_{1}$,
where $\varepsilon_{1}\sim\mathcal{N}(0,1)$, and response status after the
first stage is $X_{2,2}=I(X_{2,1}<0.7X_{1})$, which, with $A_{1}$, dictates
the feasible subset of $\mathcal{A}_{2}=\\{0,1,2,3,4\\}$ from which stage 2
treatment $A_{2}$ is assigned. Six weeks later, the outcome is generated as
$Y=\gamma_{2,0}+\gamma_{2,1}X_{1}+\gamma_{2,2}A_{1}+\gamma_{2,3}X_{2,1}+\gamma_{2,4}I(A_{2}=1)+\gamma_{2,5}I(A_{2}=2)+\gamma_{2,6}I(A_{2}=3)+\gamma_{2,7}I(A_{2}=4)+\varepsilon_{2}$,
where $\varepsilon_{2}\sim\mathcal{N}(0,1)$. With
$\boldsymbol{\gamma}_{1}=(\gamma_{1,0},\gamma_{1,1},\gamma_{1,2})^{T}=(0.00,0.90,-1.50)^{T}$
and
$\boldsymbol{\gamma}_{2}=(\gamma_{2,0},\ldots,\gamma_{2,7})=(0.00,0.30,-0.75,0.60,-0.25,-0.75,-0.75,-0.85)^{T}$,
for the $m=8$ embedded regimes defined in Figure 1,
$\\{\mathcal{V}(\mathbf{d}^{1}),\ldots,\mathcal{V}(\mathbf{d}^{8})\\}=(\theta_{1},\ldots,\theta_{8})=(-0.126,-0.374,-0.500,-0.251,-2.408,$
$-2.401,-2.494,-2.501)$, so that regime 8 is optimal, as larger reductions in
pain are desirable. The burn-in period ends at the time $t^{*}$ when each of
the $m=8$ regimes has at least 25 subjects who have completed the trial with
experience consistent with following the regime.

We compare the performance of up-front RAR using TS with $c_{t}=0.5$ and 1 for
all $t$ based on the WIPW and WAIPW estimators (1) and (2) and the unweighted
IAIPW estimator to that of simple, uniform randomization (SR) among the eight
embedded regimes. For each TS variant, at any $t$ we impose clipping constants
of 0.05 and 0.95, so if $r_{t}^{j}<0.05$ for any $j=1,\ldots,8$, set
$r_{t}^{j}=0.05$ (and likewise for 0.95), and then normalize $r_{t}^{j}$,
$j=1,\ldots,8$, to sum to one. In the WAIPW and IAIPW estimators,
$Q_{2}(\overline{\boldsymbol{x}}_{2},\overline{\boldsymbol{a}}_{2};\boldsymbol{\beta}_{2})=\beta_{2,0}+\beta_{2,1}x_{1}+\beta_{2,2}I(a_{1}=1)+\beta_{2,3}x_{2,1}+\beta_{2,4}I(a_{2}=1)+\beta_{2,5}I(a_{2}=2)+\beta_{2,6}I(a_{2}=3)+\beta_{2,7}I(a_{2}=4)$
and
$Q_{1}^{j}(x_{1},a_{1};\boldsymbol{\beta}_{1}^{j})=\beta_{1,0}^{j}+\beta_{1,1}^{j}x_{1}+\beta_{1,2}^{j}a_{1}+\beta_{1,3}^{j}x_{1}a_{1}$.

Results are displayed in Table 1. To evaluate how TS improves in-trial
outcomes, we report Monte Carlo average outcome across individuals in the
trial; average proportion of subjects assigned $A_{1}=1$, the optimal first-
stage treatment; and average proportion of subjects who were consistent with
the optimal regime within the trial. Outcomes are improved using TS over SR,
with TS resulting in lower average outcome and higher rates of assigning
optimal stage 1 treatment and optimal regime. More aggressive adaptation,
$c_{t}=1$, yields more favorable results than $c_{t}=0.5$ when using any
estimator. Results are more favorable for WAIPW- and IAIPW-based
randomization, as those estimators exploit covariate information, with the
gains mostly at stage 2. Using weighted vs. unweighted estimators yields
modest in-trial improvements.

Table 1: Simulation results using up-front RAR based on TS for 5000 Monte
Carlo (MC) replications for the scenario in Section 5.1. Columns indicate the
randomization method: WAIPW(0.5) is TS based on the WAIPW estimator (2) with
$c_{t}=0.5$ for all $t$ and AIPW(1) uses $c_{t}=1$ for all $t$; WIPW($c_{t}$)
and IAIPW($c_{t}$) are defined similarly. SR denotes simple, uniform
randomization. Mean $Y$ denotes the MC average mean outcome for the 1000
individuals in the trial; lower mean outcomes are more favorable. Proportion
$A_{1}$ Opt is the MC average proportion of subjects assigned the optimal
treatment at the first-stage; Proportion Regime Opt is the MC average
proportion of subjects who followed the optimal regime in the trial. For
estimation results, Regime Correct denotes the proportion of trials we
correctly estimate $\mathbf{d}^{8}$ to be the optimal regime;
$\mathcal{V}(\mathbf{d}^{8})$ MSE is the MC mean squared error for regime 8\.
For the estimation results, the term in parentheses, e.g., (IPW), denotes the
estimator used after the trial is completed. $\mathcal{V}(\mathbf{d}^{8})$
$95\%$ CI is the MC proportion of $95\%$ confidence intervals that cover the
true value; the term in parentheses is the estimator used to construct the
confidence interval. $\mathcal{V}(\mathbf{d}^{8})$ $95\%$ LB for lower
confidence bounds $\mathcal{V}(\mathbf{d}^{8})$ $95\%$ UB for upper confidence
bounds are defined similarly. Bold values indicate the most favorable result
among the randomization methods. Standard deviations of entries are in
parentheses

| SR | WIPW(0.5) | WIPW(1) | WAIPW(0.5) | WAIPW(1) | IAIPW(0.5) | IAIPW(1)
---|---|---|---|---|---|---|---
In Trial |  |  |  |  |  |  |
Mean $Y$ | -1.380 (0.001) | -1.795 (0.001) | -1.992 (0.001) | -1.794 (0.001) |  -1.996 (0.001) | -1.793 (0.001) | -1.993 (0.001)
Proportion $A_{1}$ Opt | 0.500 (0.000) | 0.691 (0.000) |  0.782 (0.000) | 0.691 (0.000) |  0.782 (0.000) | 0.690 (0.000) |  0.782 (0.000)
Proportion Regime Opt | 0.250 (0.000) | 0.390 (0.001) | 0.470 (0.001) | 0.401 (0.001) |  0.498 (0.002) | 0.402 (0.001) | 0.491 (0.001)
Estimation |  |  |  |  |  |  |
Regime Correct (IPW) | 0.433 (0.007) |  0.462 (0.007) | 0.444 (0.007) | 0.449 (0.007) | 0.409 (0.007) | 0.441 (0.007) | 0.430 (0.007)
Regime Correct (WIPW) | 0.397 (0.007) |  0.468 (0.007) | 0.460 (0.007) | 0.463 (0.007) | 0.436 (0.007) | 0.452 (0.007) | 0.460 (0.007)
Regime Correct (AIPW) | 0.529 (0.007) |  0.562 (0.007) | 0.534 (0.007) | 0.559 (0.007) | 0.523 (0.007) | 0.555 (0.007) | 0.534 (0.007)
Regime Correct (WAIPW) | 0.516 (0.007) | 0.548 (0.007) | 0.540 (0.007) |  0.555 (0.007) | 0.524 (0.007) | 0.545 (0.007) | 0.539 (0.007)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (IPW) | 0.814 (0.017) | 0.619 (0.013) | 0.735 (0.024) |  0.591 (0.013) | 0.671 (0.019) | 0.619 (0.013) | 0.669 (0.016)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (WIPW) | 1.185 (0.024) | 0.598 (0.013) | 0.646 (0.012) |  0.560 (0.012) | 0.576 (0.015) | 0.578 (0.012) | 0.574 (0.014)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (AIPW) | 0.544 (0.011) | 0.419 (0.008) | 0.463 (0.012) |  0.407 (0.008) | 0.446 (0.011) | 0.438 (0.009) | 0.432 (0.009)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (WAIPW) | 0.552 (0.012) | 0.411 (0.008) | 0.419 (0.010) |  0.400 (0.008) | 0.420 (0.011) | 0.427 (0.009) | 0.411 (0.009)
Coverage |  |  |  |  |  |  |
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (IPW) | 0.949 (0.003) | 0.950 (0.003) | 0.949 (0.003) | 0.954 (0.003) | 0.947 (0.003) | 0.951 (0.003) | 0.950 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (WIPW) | 0.943 (0.003) | 0.948 (0.003) | 0.945 (0.003) | 0.953 (0.003) | 0.944 (0.003) | 0.949 (0.003) | 0.948 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (AIPW) | 0.944 (0.003) | 0.952 (0.003) | 0.951 (0.003) | 0.952 (0.003) | 0.954 (0.003) | 0.950 (0.003) | 0.950 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (WAIPW) | 0.947 (0.003) | 0.951 (0.003) | 0.950 (0.003) | 0.951 (0.003) | 0.948 (0.003) | 0.947 (0.003) | 0.944 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (IPW) | 0.952 (0.003) | 0.960 (0.003) | 0.960 (0.003) | 0.954 (0.003) | 0.954 (0.003) | 0.950 (0.003) | 0.953 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (WIPW) | 0.946 (0.003) | 0.954 (0.003) | 0.948 (0.003) | 0.950 (0.003) | 0.947 (0.003) | 0.947 (0.003) | 0.945 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (AIPW) | 0.944 (0.003) | 0.959 (0.003) | 0.955 (0.003) | 0.951 (0.003) | 0.954 (0.003) | 0.950 (0.003) | 0.952 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (WAIPW) | 0.946 (0.003) | 0.957 (0.003) | 0.950 (0.003) | 0.947 (0.003) | 0.947 (0.003) | 0.944 (0.003) | 0.945 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (IPW) | 0.950 (0.003) | 0.944 (0.003) | 0.940 (0.003) | 0.955 (0.003) | 0.942 (0.003) | 0.949 (0.003) | 0.944 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (WIPW) | 0.946 (0.003) | 0.947 (0.003) | 0.944 (0.003) | 0.955 (0.003) | 0.946 (0.003) | 0.952 (0.003) | 0.947 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (AIPW) | 0.952 (0.003) | 0.947 (0.003) | 0.948 (0.003) | 0.954 (0.003) | 0.949 (0.003) | 0.944 (0.003) | 0.948 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (WAIPW) | 0.952 (0.003) | 0.946 (0.003) | 0.952 (0.003) | 0.956 (0.003) | 0.949 (0.003) | 0.948 (0.003) | 0.950 (0.003)

To evaluate post-trial performance, Table 1 shows the proportion of trials
where $\mathbf{d}^{8}$ is correctly identified as optimal and the mean-squared
error (MSE) for $\widehat{\theta}^{8}$ based on the final data at the end of
the trial. Relative to SR, using RAR based on any of the estimators generally
identifies optimal regimes at higher rates. Performance of 95% confidence
intervals and lower/upper confidence bounds for the true value based on the
asymptotic theory in Section 4 is presented for $\mathcal{V}(\mathbf{d}^{8})$;
that for other regimes is qualitatively similar. Regardless of randomization
scheme, the nominal level is achieved in almost every case. Figure 2 shows the
distribution of the 5000 centered and scaled final value estimates
$\widehat{\theta}^{j}$, $j=1,8$, obtained using the IPW, AIPW, WIPW, and WAIPW
estimators based on RAR using the WAIPW estimator with $c_{t}=1$; those for
other schemes are similar. As found by Zhang et al. (2021b) and others,
unweighted estimators result in mildly skewed distributions, while weighted
estimators yield approximate standard normality. Overall, the results suggest
that the weighted estimators effect the desired adjustment to yield asymptotic
normality. These observations hold over all simulation scenarios we have tried
with continuous and binary outcomes, including where several regimes have the
same or similar values, under which departures from normality may be
particularly problematic (e.g., Hadad et al., 2021).

Figure 2: Monte Carlo distributions of centered and scaled estimates as in the
theory of Section 4 for selected estimators in the simulation in Section 5.1.
The histograms correspond to the indicated estimator for
$\mathcal{V}(\mathbf{d}^{1})$ and $\mathcal{V}(\mathbf{d}^{8})$ under up-front
randomization using TS based on the WAIPW estimator (2) with $c_{t}=1$ for all
$t$. The vertical line indicates mean zero, and the density of a standard
normal distribution is superimposed.

###  5.2 Sequential randomization at each stage based on optimal regimes

We report on a simulation study involving 5000 Monte Carlo trials under the
two-stage scenario in Section 5.1 to compare the Q-learning-based sequential
RAR approach using TS in Section 3.2 with $c_{t}=0.25,0.50,0.75$ and 1 to SR
and two SMART-AR methods proposed by Cheung et al. (2014), a conservatively-
tuned (AR-1) and a more aggressive (AR-2) version. . To implement all RAR
methods, we posit linear models
$Q_{2}(\overline{\boldsymbol{x}}_{2},\overline{\boldsymbol{a}}_{2};\boldsymbol{\beta}_{2})=\beta_{2,0}+\beta_{2,1}x_{1}+\beta_{2,2}I(a_{1}=1)+\beta_{2,3}x_{2,1}+\beta_{2,4}I(a_{2}=1)+\beta_{2,5}I(a_{2}=2)+\beta_{2,6}I(a_{2}=3)+\beta_{2,7}I(a_{2}=4)$
and
$Q_{1}(x_{1},a_{1};\boldsymbol{\beta}_{1})=\beta_{1,0}+\beta_{1,1}x_{1}+\beta_{1,2}a_{1}$.
For sequential RAR, we set $B_{1}=b_{2}\times b_{1}=32\times 32=1024$ and
$B_{2}=1000$. The tuning parameters for AR-1 are $b=10,\tau=0.5$ and
$b=100,\tau=0.025$ for AR-2 and $\lambda_{t}=t^{-1}\tau^{(1-b)}$ for both; see
Cheung et al. (2014). Clipping constants of 0.05 and 0.95 were imposed on all
methods.

At each week, for each RAR method, newly-enrolled subjects are assigned stage
1 treatment using the same randomization probability. Already-enrolled
subjects who have reached stage 2 at this week and require stage 2
randomization are partitioned into four groups based on
$(a_{1},x_{2,2})=(0,0),(0,1),(1,0),(1,1)$. Within each group, randomization
probabilities are calculated; thus, second-stage probabilities are specific to
each stage 1 treatment-response status combination. For both AR methods, we
use the sample average of $X_{1}$ to calculate the components that make up the
stage 2 probabilities; see Cheung et al. (2014). The burn-in period ends at
the time $t^{*}$ when at least 25 subjects have completed the trial with
experience consistent with each of the $m=8$ embedded regimes.

Table 2 presents the results. As in Section 5.1, the sequential RAR method
results in improved (over SR) in-trial outcomes on average by assigning
optimal treatments and regimes at higher rates. The AR methods also improve on
SR but are relatively conservative. We were unable to achieve more favorable
results with differently-tuned versions of the AR method. For post-trial
estimation, less-aggressive RAR identifies the optimal regime at higher rates
than SR. The proposed method with $c_{t}=1$ and 0.75 yields the best post-
trial performance using any of the IPW, WIPW, AIPW, or WAIPW estimators but
lower rate of identifying the optimal regime. As expected, the AIPW and WAIPW
estimators are more efficient than the IPW and WIPW estimators. While the
primary goal of the weighted estimators is to attain nominal coverage, an
additional benefit is higher rates of identifying the optimal regime.

Table 2: Simulation results using sequential RAR based on TS for 5000 Monte
Carlo replications for the scenario in Section 5.1. Columns indicate the
randomization method: TS(0.25) is TS via Q-learning with $c_{t}=0.25$ for all
$t$, TS(0.50), TS(0.75), and TS(1) are defined similarly. AR-1 is the
conservatively-tuned version of SMART-AR and AR-2 is the more aggressive
version. SR denotes simple, uniform randomization. All entries are defined as
in Table 1.

| SR | TS(0.25) | TS(0.50) | TS(0.75) | TS(1) | AR-1 | AR-2
---|---|---|---|---|---|---|---
In Trial |  |  |  |  |  |  |
Mean $Y$ | -1.380 (0.001) | -1.976 (0.001) | -1.999 (0.001) | -2.014 (0.001) |  -2.206 (0.001) | -1.950 (0.001) | -1.957 (0.001)
Proportion $A_{1}$ Opt | 0.500 (0.000) | 0.772 (0.001) | 0.780 (0.001) | 0.785 (0.001) |  0.790 (0.000) | 0.775 (0.000) | 0.775 (0.000)
Proportion Regime Opt | 0.250 (0.000) | 0.445 (0.001) | 0.491 (0.001) | 0.517 (0.001) |  0.538 (0.002) | 0.320 (0.001) | 0.351 (0.001)
Estimation |  |  |  |  |  |  |
Regime Correct (IPW) | 0.433 (0.007) | 0.471 (0.007) | 0.463 (0.007) | 0.423 (0.007) | 0.404 (0.007) | 0.468 (0.007) |  0.484 (0.007)
Regime Correct (WIPW) | 0.397 (0.007) | 0.480 (0.007) | 0.471 (0.007) | 0.452 (0.007) | 0.434 (0.007) | 0.480 (0.007) |  0.488 (0.007)
Regime Correct (AIPW) | 0.529 (0.007) | 0.568 (0.007) | 0.549 (0.007) | 0.512 (0.007) | 0.484 (0.007) | 0.546 (0.007) |  0.569 (0.007)
Regime Correct (WAIPW) | 0.516 (0.007) | 0.569 (0.007) | 0.548 (0.007) | 0.530 (0.007) | 0.519 (0.007) | 0.558 (0.007) |  0.582 (0.007)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (IPW) | 0.814 (0.017) | 0.556 (0.012) |  0.541 (0.012) | 0.586 (0.015) | 0.660 (0.020) | 0.837 (0.018) | 0.676 (0.014)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (WIPW) | 1.185 (0.024) | 0.504 (0.011) |  0.480 (0.011) | 0.515 (0.014) | 0.517 (0.014) | 0.767 (0.016) | 0.636 (0.013)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (AIPW) | 0.544 (0.011) | 0.396 (0.009) |  0.383 (0.008) | 0.401 (0.010) | 0.433 (0.013) | 0.521 (0.011) | 0.442 (0.009)
$\mathcal{V}(\mathbf{d}^{8})$ MSE $\times 10^{2}$ (WAIPW) | 0.552 (0.012) | 0.377 (0.008) |  0.368 (0.008) | 0.381 (0.010) | 0.378 (0.010) | 0.501 (0.010) | 0.432 (0.009)
Coverage |  |  |  |  |  |  |
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (IPW) | 0.949 (0.003) | 0.940 (0.003) | 0.951 (0.003) | 0.953 (0.003) | 0.954 (0.003) | 0.949 (0.003) | 0.948 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (WIPW) | 0.943 (0.003) | 0.940 (0.003) | 0.951 (0.003) | 0.952 (0.003) | 0.954 (0.003) | 0.949 (0.003) | 0.947 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (AIPW) | 0.944 (0.003) | 0.947 (0.003) | 0.950 (0.003) | 0.958 (0.003) | 0.954 (0.003) | 0.950 (0.003) | 0.952 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% CI (WAIPW) | 0.947 (0.003) | 0.942 (0.003) | 0.947 (0.003) | 0.951 (0.003) | 0.952 (0.003) | 0.948 (0.003) | 0.951 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (IPW) | 0.952 (0.003) | 0.948 (0.003) | 0.952 (0.003) | 0.955 (0.003) | 0.960 (0.003) | 0.940 (0.003) | 0.948 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (WIPW) | 0.946 (0.003) | 0.946 (0.003) | 0.948 (0.003) | 0.950 (0.003) | 0.954 (0.003) | 0.946 (0.003) | 0.949 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (AIPW) | 0.944 (0.003) | 0.948 (0.003) | 0.948 (0.003) | 0.955 (0.003) | 0.957 (0.003) | 0.944 (0.003) | 0.949 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% LB (WAIPW) | 0.946 (0.003) | 0.946 (0.003) | 0.945 (0.003) | 0.947 (0.003) | 0.950 (0.003) | 0.946 (0.003) | 0.951 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (IPW) | 0.950 (0.003) | 0.945 (0.003) | 0.950 (0.003) | 0.950 (0.003) | 0.946 (0.003) | 0.952 (0.003) | 0.950 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (WIPW) | 0.946 (0.003) | 0.947 (0.003) | 0.951 (0.003) | 0.951 (0.003) | 0.952 (0.003) | 0.949 (0.003) | 0.948 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (AIPW) | 0.952 (0.003) | 0.946 (0.003) | 0.951 (0.003) | 0.952 (0.003) | 0.952 (0.003) | 0.954 (0.003) | 0.951 (0.003)
$\mathcal{V}(\mathbf{d}^{8})$ 95% UB (WAIPW) | 0.952 (0.003) | 0.949 (0.003) | 0.953 (0.003) | 0.957 (0.003) | 0.956 (0.003) | 0.952 (0.003) | 0.947 (0.003)

Figure 3 shows the distributions of the 5000 centered and scaled final value
estimates $\widehat{\theta}^{j}$, $j=1,8$, based on sequential RAR using TS
with $c_{t}=1$ (plots for other methods and regimes are similar). For
$\mathcal{V}(\mathbf{d}^{8})$, the weighted estimators are approximately
normal while the unweighted estimators are slightly left skewed; for
$\mathcal{V}(\mathbf{d}^{1})$, the unweighted estimators are nonnormal, and
the weighted estimators are improved, although the variance is larger than
one. We attribute this behavior to undersampling of $\mathbf{d}^{1}$, the
least effective regime; the issue is less pronounced for less aggressive
randomization. Coverage of confidence intervals and bounds for
$\mathcal{V}(\mathbf{d}^{8})$ is mildly improved for the weighted estimators.
Because of the left skew, there is notable upper confidence bound
undercoverage and lower bound overcoverage for the unweighted estimators.

Figure 3: Monte Carlo distributions of centered and scaled estimates as in the
theory of Section 4 for selected estimators in the simulation in Section 5.2.
The histograms correspond to the indicated estimator for
$\mathcal{V}(\mathbf{d}^{1})$ and $\mathcal{V}(\mathbf{d}^{8})$ under
sequential randomization using TS with $c_{t}=1$ for all $t$. The vertical
line indicates mean zero, and the density of a standard normal distribution is
superimposed.

##  6 Discussion

We have proposed methods for RAR in SMARTs using TS, where randomization can
be up-front to regimes embedded in the trial or performed sequentially at each
stage. Simulation studies demonstrate the benefits over nonadaptive
randomization: improved outcomes for subjects in the trial, improved ability
to identify an optimal regime, and little or no effect on post-trial inference
on embedded regimes. Choice of damping constant can dramatically change the
aggressiveness of TS; thus, the specific features and goals of a SMART should
be considered when specifying this parameter. When randomization is up-front,
the choice of estimator on which to base TS has some effect on performance;
WAIPW and AIPW estimators are more aggressive than WIPW and IPW estimators.
SMART-AR methods also yield good in- and post-trial performance; however, the
tuning parameters are less intuitive and effective. For any SMART, we
recommend simulating multiple different adaptive randomization methods to see
which best aligns with the trial goals.

The weighted versions of the IPW and AIPW estimators are preferred over the
unweighted estimators for post-trial inference. Normalized weighted estimators
have sampling distribution closer to standard normal, and yield improved
coverage of confidence intervals and bounds and ability to identify optimal
embedded regimes. When baseline and intermediate subject variables that are
correlated with outcome are available, we recommend the WAIPW estimator for
post-trial inference when using RAR of any type.

## References

  * Agrawal and Goyal (2012) Agrawal, S. and Goyal, N. (2012).  Analysis of Thompson sampling for the multi-armed bandit problem.  In JMLR Workshop and Conference Proceedings, volume 23, pages 39.1–39.26. 25th Annual Conference on Learning Theory.
  * Atkinson and Biswas (2013) Atkinson, A. C. and Biswas, A. (2013).  Randomised Response-Adaptive Designs in Clinical Trials.  Chapman and Hall/CRC Press, Boca Raton, Florida.
  * Berry (2011) Berry, D. (2011).  Adaptive clinical trials in oncology.  Nature Reviews Clinical Oncology 9, 199–207.
  * Berry (2015) Berry, D. A. (2015).  The brave new world of clinical cancer research: Adaptive biomarker-driven trials integrating clinical practice with clinical research.  Molecular Oncology 9, 951–959.
  * Berry and Fristedt (1985) Berry, D. A. and Fristedt, B. (1985).  Bandit Problems: Sequential Allocation of Experiments.  Chapman and Hall, London.
  * Berry et al. (2010) Berry, S. M., Carlin, B. P., Lee, J. J., and Muller, P. (2010).  Bayesian Adaptive Methods for Clinical Trials.  Chapman and Hall/CRC Press, Boca Raton, Florida.
  * Bibaut et al. (2021) Bibaut, A., Dimakopoulou, M., Kallus, N., Chambaz, A., and van der Laan, M. (2021).  Post-contextual-bandit inference.  Advances in Neural Information Processing Systems 34, 28548–28559.
  * Bigirumurame et al. (2022) Bigirumurame, T., Uwimpuhwe, G., and Wason, J. (2022).  Sequential multiple assignment randomized trial studies should report all key components: a systematic review.  Journal of Clinical Epidemiology 142, 152–160.
  * Cheung et al. (2014) Cheung, Y. K., Chakraborty, B., and Davidson, K. W. (2014).  Sequential multiple assignment randomized trial (SMART) with adaptive randomization for quality improvement in depression treatment program.  Biometrics 71, 450–459.
  * Food and Drug Administration (2019) Food and Drug Administration (2019).  Adaptive Designs for Clinical Trials of Drugs and Biologics.  Guidance for industry, November 2019, https://www.fda.gov/media/78495/download.
  * Godambe (1985) Godambe, V. P. (1985).  The foundations of finite sample estimation in stochastic processes.  Biometrika 72, 419–428.
  * Godambe and Heyde (2010) Godambe, V. P. and Heyde, C. C. (2010).  Quasi-likelihood and optimal estimation.  International Statistical Review 55, 386–399.
  * Hadad et al. (2021) Hadad, V., Hirshberg, D. A., Zhan, R., Wager, S., and Athey, S. (2021).  Confidence intervals for policy evaluation in adaptive experiments.  Proceedings of the National Academy of Sciences 118, e2014602118.
  * Hall and Heyde (1980) Hall, P. and Heyde, C. C. (1980).  Martingale Limit Theory and its Application.  Academic Press, New York.
  * Heyde (1997) Heyde, C. C. (1997).  Quasi-Likelihood and Its Application: A General Approach to Optimal Parameter Estimation.  Springer, New York.
  * Hu and Rosenberger (2006) Hu, F. and Rosenberger, W. F. (2006).  The Theory of Response-Adaptive Randomization in Clinical Trials.  John Wiley and Sons, New York.
  * Kidwell (2014) Kidwell, K. M. (2014).  SMART designs in cancer research: Past, present, and future.  Clinical Trials 11, 445–456.
  * Kim et al. (2011) Kim, E. S., Herbst, R. S., Wistuba, I. I., Lee, J. J., Blumenschein, G. R., Tsao, A., Stewart, D. J., Hicks, M. E., Erasmus, J., Gupta, S., et al. (2011).  The BATTLE trial: personalizing therapy for lung cancer.  Cancer Discovery 1, 44–53.
  * Laber et al. (2014) Laber, E. B., Lizotte, D. J., Qian, M., Pelham, W. E., and Murphy, S. A. (2014).  Dynamic treatment regimes: Technical challenges and applications.  Electronic Journal of Statistics 8,.
  * Lattimore and Szepesvári (2020) Lattimore, T. and Szepesvári, C. (2020).  Bandit Algorithms.  Cambridge University Press.
  * Lavori and Dawson (2004) Lavori, P. W. and Dawson, R. (2004).  Dynamic treatment regimes: Practical design considerations.  Clinical Trials 1, 9–20.
  * Liu et al. (2023) Liu, X., Deliu, N., and Chakraborty, B. (2023).  Microrandomized trials: Developing just-in-time adaptive interventions for better public health.  American Journal of Public Health 113, 60–69.
  * Lorenzoni et al. (2023) Lorenzoni, G., Petracci, E., Scarpi, E., Baldi, I., Gregori, D., and Nanni, O. (2023).  Use of sequential multiple assignment randomized trials (smarts) in oncology: systematic review of published studies.  British Journal of Cancer 128, 1177–1188.
  * Manschot et al. (2023) Manschot, C., Laber, E., and Davidian, M. (2023).  Interim monitoring of sequential multiple assignment randomized trials using partial information.  Biometrics 79, 2881–2894.
  * Moodie et al. (2014) Moodie, E. E., Dean, N., and Sun, Y. R. (2014).  Q-learning: Flexible learning about useful utilities.  Statistics in Biosciences 6, 223–243.
  * Murphy (2005) Murphy, S. A. (2005).  An experimental design for the development of adaptive treatment strategies.  Statistics in Medicine 24, 1455–1481.
  * Proschan and Evans (2020) Proschan, M. and Evans, S. (2020).  Resist the temptation of response-adaptive randomization.  Clinical Infectious Diseases 71, 3002–3004.
  * Russo et al. (2018) Russo, D., Roy, B. V., Kazerouni, A., Osband, I., and Wen, Z. (2018).  A tutorial on Thompson sampling.  Foundations and Trends in Machine Learning 11, 1–96.
  * Sheu and Yao (1984) Sheu, S. S. and Yao, Y. S. (1984).  A strong law of large numbers for martingales.  Proceedings of the American Mathematical Society 92, 283–287.
  * Somers et al. (2023) Somers, T. J., Winger, J. G., Fisher, H. M., Hyland, K. A., Davidian, M., Laber, E. B., and et al. (2023).  Behavioral cancer pain intervention dosing: results of a sequential multiple assignment randomized trial.  PAIN 164, 1935–1941.
  * Tewari and Murphy (2017) Tewari, A. and Murphy, S. A. (2017).  From ads to interventions: Contextual bandits in mobile health.  In Rehg, J. M., Murphy, S. A., and Kuman, S., editors, Mobile Health: Sensors, Analytic Methods, and Applications, pages 495–517, New York. Springer.
  * Thall and Wathen (2007) Thall, P. and Wathen, K. (2007).  Practical Bayesian adaptive randomization in clinical trials.  European Journal of Cancer 43, 859–866.
  * Thall et al. (2015) Thall, P. F., Fox, P. S., and Wathen, J. K. (2015).  Some caveats for outcome adaptive randomization in clinical trials.  In Sverdiov, O., editor, Modern Adaptive Randomized Clinical Trials: Statistical and Practical Aspects, pages 287–305, Boca Raton, Florida. Chapman and Hall/CRC Press.
  * Thompson (1933) Thompson, W. R. (1933).  On the likelihood that one unknown probability exceeds another in view of the evidence of two samples.  Biometrika 25, 285–294.
  * Tsiatis (2006) Tsiatis, A. A. (2006).  Semiparametric Theory and Missing Data.  Springer, New York.
  * Tsiatis et al. (2020) Tsiatis, A. A., Davidian, M., Holloway, S., and Laber, E. B. (2020).  Dynamic Treatment Regimes: Statistical Methods for Precision Medicine.  Chapman and Hall/CRC Press, Boca Raton, Florida.
  * Viele et al. (2020) Viele, K., Broglio, K., McGlothlin, A., and Saville, B. R. (2020).  Comparison of methods for control allocation in multiple arm studies using response adaptive randomization.  Clinical Trials 17, 52–60.
  * Villar et al. (2015) Villar, S. S., Bowden, J., and Wason, J. (2015).  Multi-armed bandit models for the optimal design of clinical trials: Benefits and challenges.  Statistical Science 30, 199–215.
  * Villar et al. (2021) Villar, S. S., Robertson, D. S., and Rosenberger, W. F. (2021).  The temptation of overgeneralizing response-adaptive randomization.  Clinical Infectious Diseases 73, e842.
  * Wang and Yee (2019) Wang, H. and Yee, D. (2019).  I-SPY 2: a neoaduvant adaptive clinical trial designed to improve outcomes in high-risk breast cancer.  Current Breast Cancer Reports 11, 303–310.
  * Wang et al. (2022) Wang, J., Wu, L., and Wahed, A. S. (2022).  Adaptive randomization in a two-stage sequential multiple assignment randomized trial.  Biostatistics 23, 1182–1199.
  * Xie and Singh (2013) Xie, M. and Singh, K. (2013).  Confidence distribution, the frequentist distribution estimator of a parameter: A review.  International Statistical Review 81, 3–39.
  * Zhang et al. (2021a) Zhang, K. W., Janson, L., and Murphy, S. A. (2021a).  Inference for batched bandits.  arXiv preprint arXiv:2002.03217 .
  * Zhang et al. (2021b) Zhang, K. W., Janson, L., and Murphy, S. A. (2021b).  Statistical inference with M-estimators on adaptively collected data.  Advances in Neural Information Processing Systems 34, 7460–7471.

## Appendix A: Interim AIPW Estimator

A limitation of using the WIPW or WAIPW estimators as the basis for updating
randomization probabilities is that only the data on subjects in
$\mathfrak{D}_{t-1}$ who have completed the trial and have $Y$ observed are
used. Thus, partial information from subjects still progressing through the
trial is ignored. To make use of this information, we adopt the Interim AIPW
(IAIPW) estimator of Manschot, Laber, and Davidian (2023), which is derived by
viewing the partial information through the lens of monotone coarsening
(Tsiatis, 2006).

We present the estimator in the case where $K=2$ and
$n_{t-1}=\sum^{N}_{i=1}\Gamma_{t-1,i}$ ; a general formulation is presented in
Manschot et al. (2023). In this setting, the weighted IAIPW (WIAIPW) estimator
for $\theta^{j}=\mathcal{V}(\mathbf{d}^{j})$ based on $\mathfrak{D}_{t-1}$ is

$\displaystyle\widehat{\theta}_{t}^{WIAIPW,j}=n_{t-1}^{-1}\sum^{N}_{i=1}$
$\displaystyle\Gamma_{t-1,i}W_{\tau_{i}}^{j}\left[\frac{\Delta_{t-1,i}C_{i}^{j}Y_{i}}{\pi_{\tau_{i},2}^{j}(\overline{\boldsymbol{X}}_{2,i})\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})\widehat{\nu}_{\Delta,t-1}}+\left\\{1-\frac{\overline{C}_{1,i}^{j}I(\kappa_{t-1,i}=2)}{\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})\widehat{\nu}_{2,t-1}}\right\\}L_{1}^{j}(\boldsymbol{X}_{1,i})\right.$
$\displaystyle+\left.\left\\{\frac{\overline{C}_{1,i}^{j}I(\kappa_{t-1,i}=2)}{\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})\widehat{\nu}_{2,t-1}}-\frac{\Delta_{t-1,i}C_{i}^{j}}{\pi_{\tau_{i},2}^{j}(\overline{\boldsymbol{X}}_{2,i})\pi_{\tau_{i},1}^{j}(\boldsymbol{X}_{1,i})\widehat{\nu}_{\Delta,t-1}}\right\\}L^{j}_{2}(\overline{\boldsymbol{X}}_{2,i})\right],$

where $\widehat{\nu}_{\Delta,t-1}=\sum^{N}_{i=1}\Delta_{t-1,i}/n_{t-1}$ and
$\widehat{\nu}_{2,t-1}=\sum^{N}_{i=1}\Gamma_{t-1,i}I(\kappa_{t-1,i}=2)/n_{t-1}$
are estimators for $P(\Delta_{t-1}=1|\Gamma_{t-1}=1,\mathfrak{D}_{t-1})$ and
$P(\kappa_{t-1}=2|\Gamma_{t-1}=1,\mathfrak{D}_{t-1})$, $W_{\tau_{i}}^{j}$ is a
weight, and $L^{j}_{k}(\overline{\boldsymbol{x}}_{k})$, $k=1,2$, are arbitrary
functions of $\overline{\boldsymbol{x}}_{k}$. The efficient choice for the
functions $L^{j}_{k}(\overline{\boldsymbol{x}}_{k}),\,k=1,2$ are derived from
the $Q$-functions which can be modeled and estimated based on a backward
iterative scheme similar to that presented for the WAIPW estimator, with
modifications to the definition of the pseudo outcomes for subjects in
$\mathfrak{D}_{t-1}$ who have not yet completed the trial; see Manschot et al.
(2023). As for (1) and (2), the sampling distribution of
$\widehat{\boldsymbol{\theta}}^{WIAIPW}_{t}=(\widehat{\theta}_{t}^{WIAIPW,1},\ldots,\widehat{\theta}_{t}^{WIAIPW,m})^{T}$
can be approximated using large sample theory analogous to that presented in
the following two sections. The general form of the weights is complex and in
preliminary simulations, weighted versions of the estimators did not
outperform a version where the weights are one for all $\tau$. Due to this,
the simulations in Section 5 only use the unweighted IAIPW as a basis for
randomization for simplicity. Moreover, we did not find the added precision
from using partial information to have a notable impact in randomization.

## Appendix B: Theory on Weighted Estimators

Overview

We use the theory for Martingale estimating functions (Godambe, 1985; Heyde,
1997; Godambe and Heyde, 2010, MEFs) to characterize the asymptotic behavior
or our post-trial estimator under Thompson Sampling. For completeness, we
state and prove the results we will use.

Suppose that we observe the first $T$ elements of a discrete-time stochastic
process, $\bm{Z}_{t},\ldots,\bm{Z}_{T}$, taking values in $\mathcal{Z}$. In
addition, suppose that the law of this stochastic process is $P$ and our
interest is in estimating
$\boldsymbol{\theta}^{*}=\boldsymbol{\theta}(P)\in\bm{\Theta}\subseteq\mathbb{R}^{p}$.
We assume there exists an $\mathbb{R}^{p}$-valued process on $\bm{\Theta}$
given by

$\mathcal{M}_{T}(\boldsymbol{\theta})=\sum_{t=1}^{T}\boldsymbol{M}_{t}(\boldsymbol{\theta}),$

such that $\left\\{\boldsymbol{M}_{t}(\boldsymbol{\theta}^{*})\right\\}_{t\geq
1}$ is a Martingale difference sequence with respect to the filtration
$\left\\{\mathcal{F}_{t}\right\\}_{t\geq 1}$, i.e.,
$E\left\\{\boldsymbol{M}_{t}(\boldsymbol{\theta}^{*})\big{|}\mathcal{F}_{t-1}\right\\}=0$
with probability one for all $t\geq 1$. We construct an estimator
$\widehat{\boldsymbol{\theta}}_{T}$ as a root of
$\mathcal{M}_{T}(\boldsymbol{\theta})=0$; more generally, it could be that no
exact root exists for finite $T$ and thus we select
$\widehat{\boldsymbol{\theta}}_{T}$ as a minimizer of
$||\mathcal{M}_{T}(\boldsymbol{\theta})||_{2}$, where $||\cdot||_{2}$ denotes
the Euclidean norm.

Consistency

We first consider conditions under which
$\widehat{\boldsymbol{\theta}}_{T}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\boldsymbol{\theta}^{*}$. Define
$\bm{\xi}_{T}(\boldsymbol{\theta})=\sum_{t=1}^{T}\mathrm{var}\left\\{\boldsymbol{M}_{t}(\boldsymbol{\theta})\big{|}\mathcal{F}_{t-1}\right\\}$
We will make use of the following conditions.

  * (A0)

$\left\\{\boldsymbol{M}_{t}(\boldsymbol{\theta}^{*})\right\\}_{t\geq 1}$ is a
Martingale difference sequence with respect to
$\left\\{\mathcal{F}_{t}\right\\}_{t\geq 1}$.

  * (A1)

$\mathcal{M}_{T}(\boldsymbol{\theta})$ is continuously differentiable for all
$\boldsymbol{\theta}\in\bm{\Theta}$ and
$\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\boldsymbol{\theta})$ is
invertible with probability one for all $\boldsymbol{\theta}\in\bm{\Theta}$
provided $T$ is sufficiently large.

  * (A2)

The minimum eigenvalue of $\bm{\xi}_{T}(\boldsymbol{\theta}^{*})$ satisfies
$\lambda_{\min}\left\\{\bm{\xi}_{T}(\boldsymbol{\theta}^{*})\right\\}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\infty,$ as $T\rightarrow\infty$.

  * (A3)

For any $\boldsymbol{\theta}\in\bm{\Theta}$

$\sup_{\boldsymbol{\theta}\in\bm{\Theta}}\bigg{|}\bigg{|}\left\\{\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{t}(\boldsymbol{\theta})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}\bigg{|}\bigg{|}=O_{P}(1).$

  * (A4)

$\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\boldsymbol{\theta}^{*})=O_{P}(1)$.

###### Theorem 1

Assume (A0)-(A4) and suppose that $\widehat{\boldsymbol{\theta}}_{t}$
satisfies $\mathcal{M}_{t}(\widehat{\boldsymbol{\theta}}_{t})=o_{P}(1)$, then
$\widehat{\boldsymbol{\theta}}_{T}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\boldsymbol{\theta}^{*}$ as $T\rightarrow\infty$.

###### Proof.

Using a Taylor’s expansion write

$\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})=\mathcal{M}_{T}(\boldsymbol{\theta}^{*})+\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\left(\widehat{\boldsymbol{\theta}}_{T}-\boldsymbol{\theta}^{*}\right),$

where $\overline{\boldsymbol{\theta}}_{T}$ is an intermediate point between
$\widehat{\boldsymbol{\theta}}_{T}$ and $\boldsymbol{\theta}^{*}$. Using (A1),
we can re-arrange the above expression to obtain

$\displaystyle\widehat{\boldsymbol{\theta}}_{T}$ $\displaystyle=$
$\displaystyle\boldsymbol{\theta}^{*}+\left\\{\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\right\\}^{-1}\left\\{\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})-\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\right\\}$
$\displaystyle=$
$\displaystyle\boldsymbol{\theta}^{*}+\left\\{\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}\left\\{\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})-\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\right\\}$
$\displaystyle=$
$\displaystyle\boldsymbol{\theta}^{*}+\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\left\\{\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\left\\{\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})-\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\right\\}.$

From the minimum eigenvalue condition
$||\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})||_{F}=o_{P}(1)$, where
$||\cdot||_{F}$ is the Frobenius norm. Furthermore, from (A3) it follows that

$\big{|}\big{|}\left\\{\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}\big{|}\big{|}=O_{P}(1).$

Finally, we see that
$\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})$
is the product of $o_{P}(1)$ terms and thus is $o_{p}(1)$, and
$\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\boldsymbol{\theta}^{*})=O_{P}(1)$
from (A4). Thus, we have shown
$\widehat{\boldsymbol{\theta}}_{T}=\boldsymbol{\theta}^{*}+o_{P}(1)$. ∎

Asymptotic Normality

To establish normality we replace (A4) with

  * (A5)

Let $\boldsymbol{I}_{p}$ denote the $p\times p$ identity matrix,
$\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\stackrel{{\scriptstyle
D}}{{\longrightarrow}}\mathcal{N}(0,\boldsymbol{I}_{p})$.

In addition, we assume the following conditions.

  * (A6)

For any consistent estimator
$\widetilde{\boldsymbol{\theta}}_{T}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\boldsymbol{\theta}^{*}$

$\bigg{|}\bigg{|}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\left\\{\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\widetilde{\boldsymbol{\theta}}_{T})-\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\right\\}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\bigg{|}\bigg{|}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0,$

as $T\rightarrow\infty$.

  * (A7)

There exists fixed and positive definite matrix $\bm{\Sigma}$ such that

$-\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\bm{\Sigma},$

as $T\rightarrow\infty$.

  * (A8)

$\widehat{\boldsymbol{\theta}}_{T}$ satisfies
$\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0$ as $T\rightarrow\infty$.

###### Theorem 2

Assume (A0)-(A3) and(A5)-(A8), then
$\bm{\Sigma}\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})\left(\widehat{\boldsymbol{\theta}}_{T}-\boldsymbol{\theta}^{*}\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}\mathcal{N}(0,\boldsymbol{I}_{p})$ as
$T\rightarrow\infty$.

###### Proof.

As in our proof of consistency, we begin with a Taylor series expansion

$\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})=\mathcal{M}_{T}(\boldsymbol{\theta}^{*})+\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\left(\widehat{\boldsymbol{\theta}}_{T}-\boldsymbol{\theta}^{*}\right),$

where $\overline{\boldsymbol{\theta}}_{T}$ is an intermediate point between
$\widehat{\boldsymbol{\theta}}_{T}$ and $\boldsymbol{\theta}^{*}$. re-
arranging terms yields

$\displaystyle(\widehat{\boldsymbol{\theta}}_{T}-\boldsymbol{\theta}^{*})$
$\displaystyle=$
$\displaystyle\left\\{\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\right\\}^{-1}\left\\{\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})-\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\right\\}$
$\displaystyle\Leftrightarrow\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})(\widehat{\boldsymbol{\theta}}_{T}-\boldsymbol{\theta}^{*})$
$\displaystyle=$
$\displaystyle\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})\left\\{\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}$
$\displaystyle\times$
$\displaystyle\left\\{\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})-\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\right\\}$
$\displaystyle=$
$\displaystyle\left\\{\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\left\\{\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})-\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\right\\}$
$\displaystyle=$
$\displaystyle\left\\{-\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\boldsymbol{\theta}^{*})$
$\displaystyle\quad+\left\\{\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})$
$\displaystyle=$ $\displaystyle\bm{Z}_{T}+\boldsymbol{E}_{T}.$

We show that $\bm{Z}_{T}\stackrel{{\scriptstyle
D}}{{\longrightarrow}}\mathcal{N}(0,\bm{\Sigma}^{-1}\bm{\Sigma}^{-T})$ and
that $\boldsymbol{E}_{T}=o_{P}(1)$, where $\bm{\Sigma}^{-T}$ is the transpose
of $\bm{\Sigma}^{-1}$. To obtain the limit for $\bm{Z}_{T}$ we write

$\displaystyle\bm{Z}_{T}$
$\displaystyle=\left\\{-\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\boldsymbol{\theta}^{*})$
$\displaystyle=\left[-\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\boldsymbol{\theta}^{*})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})+\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\left\\{\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\boldsymbol{\theta}^{*})-\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\right\\}\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right]^{-1}$
$\displaystyle\times\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\boldsymbol{\theta}^{*}),$

applying (A6), it follows that the term inside the square brackets is
$\bm{\Sigma}^{-1}+o_{P}(1)$ so that

$\bm{Z}_{T}=\bm{\Sigma}^{-1}\bm{\xi}_{T}^{1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\boldsymbol{\theta}^{*})+o_{P}(1)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}\mathcal{N}(0,\bm{\Sigma}^{-1}\bm{\Sigma}^{-T}).$

It remains to show that $\boldsymbol{E}_{T}=o_{P}(1)$. From (A3)
$\left\\{\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\nabla_{\boldsymbol{\theta}}\mathcal{M}_{T}(\overline{\boldsymbol{\theta}}_{T})\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\right\\}^{-1}=O_{P}(1)$.
Furthermore,
$\bm{\xi}_{T}^{-1/2}(\boldsymbol{\theta}^{*})\mathcal{M}_{T}(\widehat{\boldsymbol{\theta}}_{T})=o_{P}(1)$
from (A8). This proves the result. ∎

### Choosing the Weights

WIPW Estimator After Up-Front Randomization

Consider the WIPW estimator for the value $\theta^{j}$ of $j$th embedded
regime in (1) at the end of the trial under the conditions of Section 4, so
that subjects are indexed by $t$. We require that the weights $W^{j}_{t}$ be
chosen as functions of $\mathfrak{D}_{t-1}$ so that (i) the estimating
equations are conditionally unbiased,
$E\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})|\mathfrak{D}_{t-1}\\}=0$,
and (ii) the variance is stabilized,
$E[\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})\\}^{2}|\mathfrak{D}_{t-1}]$
$=\sigma>0$ for all $t$, where
$M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})$
is the estimating function (4) given by

$M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})=\frac{W^{j}_{t}C^{j}_{t}}{\\{\prod_{k=2}^{K}\pi_{t,k}^{j}(\overline{\boldsymbol{X}}_{k,t})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1,t})}(Y_{t}-\theta^{j}).$

As in Section 4, take
$\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})=\\{\boldsymbol{X}_{1},\boldsymbol{X}_{2}^{*}(\mathbf{d}^{j}),\ldots,\boldsymbol{X}_{k}^{*}(\mathbf{d}^{j})\\}$,
$k=1,\ldots,K$.

To show (i), when $C^{j}_{t}=1$, the consistency assumption given in Section
2.1 implies that $Y_{t}=Y^{*}(\mathbf{d}^{j})$ and
$\overline{\boldsymbol{X}}_{K,t}=\overline{\boldsymbol{X}}_{K}^{*}(\mathbf{d}^{j})$.
Thus, if $W^{j}_{t}$ is a function of $\mathfrak{D}_{t-1}$,

$\displaystyle E$
$\displaystyle\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})|\mathfrak{D}_{t-1}\\}$
$\displaystyle=W^{j}_{t}E\left[\left.\frac{C^{j}_{t}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}}{\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}\right|\mathfrak{D}_{t-1}\right]$
$\displaystyle=W^{j}_{t}E\left(\left.E\left[\left.\frac{C^{j}_{t}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}}{\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}\right|\mathcal{W}^{*},\mathfrak{D}_{t-1}\right]\right|\mathfrak{D}_{t-1}\right)$
$\displaystyle=W^{j}_{t}E\left[\left.\frac{E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})}{\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}\right|\mathfrak{D}_{t-1}\right]$
$\displaystyle=W^{j}_{t}E\left[\left.\frac{\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}{\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}\right|\mathfrak{D}_{t-1}\right]$
$\displaystyle=W^{j}_{t}E\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}=0,$

where we have used the result analogous to that in Section 6.4.3 of Tsiatis et
al. (2020) that, under the consistency, positivity, and sequential
ignorability assumptions,
$E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})$
and the fact that $\mathcal{W}^{*}$ and thus $Y^{*}(\mathbf{d}^{j})$ is
independent of $\mathfrak{D}_{t-1}$.

To find $W^{j}_{t}$ satisfying (ii), we have, using the consistency
assumption, $(C^{j}_{t})^{2}=C^{j}_{t}$, and that
$E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})$,
that

$\displaystyle E$
$\displaystyle\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})^{2}|\mathfrak{D}_{t-1}\\}=(W^{j}_{t})^{2}E\left(\left.\frac{C^{j}_{t}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}}{\big{[}\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})\big{]}^{2}}\right|\mathfrak{D}_{t-1}\right)$
$\displaystyle=(W^{j}_{t})^{2}E\left(\left.\frac{E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})}{\big{[}\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})\big{]}^{2}}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}\right|\mathfrak{D}_{t-1}\right)$
$\displaystyle=(W^{j}_{t})^{2}E\left[\left.\frac{\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}}{\prod_{k=2}^{K}\pi_{t,k}^{j}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}\right|\mathfrak{D}_{t-1}\right]$
(B.1)

as in 6.

To ensure (ii), we thus want to choose $W^{j}_{t}$ so that (B.1) is a constant
that can be $j$-specific but that does not depend on $\mathfrak{D}_{t-1}$. For
definiteness and simplicity, consider $K=2$; the following argument extends to
general $K$. As in Section 4 for the cancer pain management SMART,
randomization at stage 1 ordinarily does not depend on covariate information,
which implies that $\pi_{t,1}^{j}(\boldsymbol{X}_{1})=\pi^{j}_{t,1}$ depending
on $\mathfrak{D}_{t-1}$ and $j$ but not $X_{1}$. Randomization at stages
beyond stage 1 is typically to sets of feasible treatment options within each
level of a binary variable such as response status, and depends on no other
covariate information. Let $R^{*}(\mathbf{d}^{j})$ be the potential binary
response status variable taking on values 0 and 1 for a subject whose first
stage treatment assignment is consistent with regime $\mathbf{d}^{j}$, where
$R^{*}(\mathbf{d}^{j})$ is a function of or component of
$\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})$
($R^{*}(\mathbf{d}^{j})=X_{22}^{*}(\mathbf{d}^{j})$ in the cancer pain
management SMART). By the consistency assumption, the observed response status
$R$ is equal to $R^{*}(\mathbf{d}^{j})$ when $C^{j}_{t}=1$. Then, as in the
cancer pain management SMART, it follows that
$\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}=\pi^{j(s)}_{t,2}$
when $R^{*}(\mathbf{d}^{j})=s$, $s=0,1$, where for each $r$,
$\pi^{j(s)}_{t,2}$ depends on $\mathfrak{D}_{t-1}$ and $j$ but not on
$\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})$, so that

$\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}=I\\{R^{*}(\mathbf{d}^{j})=0\\}\pi^{j(0)}_{t,2}+I\\{R^{*}(\mathbf{d}^{j})=1\\}\pi^{j(1)}_{t,2}.$
(B.2)

Because $\mathcal{W}^{*}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathfrak{D}_{t-1}$, generalizing (7), it
follows that (B.1) is equal to

$(W^{j}_{t})^{2}\left(\frac{\mu^{j(1)}}{\pi^{j(1)}_{t,2}\pi_{t,1}^{j}}+\frac{\mu^{j(0)}}{\pi^{j(0)}_{t,2}\pi_{t,1}^{j}}\right),$
(B.3)

where
$\mu^{j(s)}=E[I\\{R^{*}(\mathbf{d}^{j})=s\\}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}]$,
$s=0,1$. Thus, we wish to choose $W^{j}_{t}$ so that (B.3) is equal to a
constant depending only on $j$ and not on $\mathfrak{D}_{t-1}$. Denoting this
constant as $\Xi^{j}$, the weights should be chosen so that

$W^{j}_{t}=(\Xi^{j})^{1/2}\left(\frac{\mu^{j(1)}}{\pi^{j(1)}_{t,2}\pi_{t,1}^{j}}+\frac{\mu^{j(0)}}{\pi^{j(0)}_{t,2}\pi_{t,1}^{j}}\right)^{-1/2}.$
(B.4)

Of course, in practice the constant $\Xi^{j}$ is unknown. If there is a burn-
in period of length $t^{*}$ weeks during which nonadaptive randomization takes
place with $W^{j}_{t}\equiv 1$, $t\leq t^{*}$, from (B.3), it must be the case
that $\Xi^{j}$ satisfies

$\Xi^{j}=\frac{\mu^{j(1)}}{\pi^{j(1)}_{t^{*},2}\pi_{t^{*},1}^{j}}+\frac{\mu^{j(0)}}{\pi^{j(0)}_{t^{*},2}\pi_{t^{*},1}^{j}},$
(B.5)

where the randomization probabilities during the burn-in period are fixed, so
that $\pi^{j(s)}_{t^{*},2}$, $s=0,1,$ and $\pi_{t^{*},1}^{j}$ are known
constants not depending on the data. Thus, if $\mu^{j(s)}$, $s=0,1$, were
known, so that $\Xi^{j}$ is known, the weights for $t\geq t^{*}+1$ should be
chosen as in (B.4) with $\Xi^{j}$ as in (B.5).

These observations suggest our proposed scheme for implementation in practice.
We propose to obtain weights in this spirit by estimating $\Xi^{j}$ by an
estimator $\widehat{\Xi}_{t^{*}}^{j}$, say, based on the data from a burn-in
period of length $t^{*}$ weeks during which nonadaptive randomization takes
place with $W^{j}_{t}\equiv 1$, $t\leq t^{*}$, and then treat the estimate as
fixed and known henceforth. Specifically, let $\widehat{\mu}^{j(s)}_{t}$ be
estimators for $\mu^{j(s)}$, $s=0,1$, at any time $t$ based on
$\mathfrak{D}_{t-1}$, as in the case of the cancer pain management SMART (in
the original notation) in Section 4; with the indexing scheme here (by $t$),
$\widehat{\mu}^{j(s)}_{t}=t^{-1}\sum_{u=1}^{t}I(R_{u}=s)C^{j}_{u}(Y_{u}-\widetilde{\theta}^{j}_{t})^{2}/(\pi^{j(s)}_{u,2}\pi^{j}_{u,1})$,
$s=0,1$, where $\widetilde{\theta}^{j}_{t}$ is an estimator for $\theta^{j}$
using $\mathfrak{D}_{t-1}$ (e.g., the unweighted IPW estimator). Define

$\widehat{\Xi}^{j}_{t}=\frac{\widehat{\mu}^{j(1)}_{t}}{\pi^{j(1)}_{t,2}\pi_{t,1}^{j}}+\frac{\widehat{\mu}^{j(0)}_{t}}{\pi^{j(0)}_{t,2}\pi_{t,1}^{j}}.$

Then (B.5) implies that the estimator $\widehat{\Xi}_{t^{*},j}$ for $\Xi^{j}$
based on the burn-in data is

$\widehat{\Xi}^{j}_{t^{*}}=\frac{\widehat{\mu}^{j(1)}_{t^{*}}}{\pi^{j(1)}_{t^{*},2}\pi_{t^{*},1}^{j}}+\frac{\widehat{\mu}^{j(0)}_{t^{*}}}{\pi^{j(0)}_{t^{*},2}\pi_{t^{*},1}^{j}},$

where as above $\pi^{j(s)}_{t^{*},2}$, $s=0,1$, and $\pi_{t^{*},1}^{j}$ are
known constants not depending on the data. Treating
$\widehat{\Xi}^{j}_{t^{*}}$ as fixed and known, for $t\geq t^{*}+1$, from
(B.4), take
$W^{j}_{t}=(\widehat{\Xi}^{j}_{t^{*}}/\widehat{\Xi}^{j}_{t})^{1/2}$. Because
$\widehat{\Xi}^{j}_{t^{*}}$ depends only on
$\mathfrak{D}_{t^{*}}\subseteq\mathfrak{D}_{t-1}$ for $t\geq t^{*}+1$, and
$\widehat{\Xi}^{j}_{t}$ depends only on $\mathfrak{D}_{t-1}$, the weights
$W^{j}_{t}$ for $t\geq t^{*}+1$ depend on $\mathfrak{D}_{t-1}$, as required.

WIPW Estimator After Sequential Randomization

Again consider $K=2$ for definiteness. The preceding scheme for obtaining
weights assumes that up-front randomization is used, so that, returning to the
original notation, for a subject who entered the SMART at time $\tau$, the
probabilities in the denominator of (1) for that subject’s contribution are
determined at that point based on $\mathfrak{D}_{\tau-1}$. Thus, if up-front
randomization was used, the estimators for $\mu^{(s)}_{j}$, $s=0,1$, and the
corresponding weight $W^{j}_{\tau}$ for a such a subject depend on
$\pi^{j(s)}_{\tau,2}$, $s=0,1$, and $\pi_{\tau,1}^{j}$.

If the sequential randomization scheme is used, this calculation should be
modified. In obtaining the WIPW estimator for $\theta^{j}$ at the end of the
trial, for a subject who entered the SMART at time $\tau$ and then reached
stage 2 at time $\tau+v$, say, the probabilities in the denominator of (1)
should be $\pi_{\tau,1}^{j}$ and $\pi^{j(s)}_{\tau+v,2}$, $s=0,1$. Likewise,
the estimators for $\mu^{(s)}_{j}$, $s=0,1$, and the corresponding weight
should be calculated based on $\pi^{j(s)}_{\tau+v,2}$, $s=0,1$, and
$\pi_{\tau,1}^{j}$, so that the weight for this subject depends on
$\mathfrak{D}_{\tau+v-1}$ and thus can be written as $W^{j}_{\tau+v}$. When
$\tau+v<t^{*}$, $W^{j}_{\tau+v}=1$.

WAIPW Estimator After Up-Front Randomization

Analogous to the developments for the WIPW estimator for $\theta^{j}$ for
embedded regime $j$, using the indexing scheme (by $t$) in Section 4, we
identify the corresponding estimating function
$M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})$,
demonstrate that (i)
$E\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})|\mathfrak{D}_{t-1}\\}=0$,
and determine weights $W^{A,j}_{t}$ so that (ii) the variance
$E[\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})\\}^{2}|\mathfrak{D}_{t-1}]$
is a constant depending only on $j$. We demonstrate for $K=2$; the argument
extends to general $K$.

From 2, for $K=2$ it is straightforward that

$\displaystyle M_{t}^{j}$
$\displaystyle(\overline{\boldsymbol{X}}_{2,t},\overline{\boldsymbol{A}}_{2,t},Y_{t};\theta^{j})=W^{A,j}_{t}\left[\frac{C^{j}_{t}Y_{t}}{\pi^{j}_{t,2}(\overline{\boldsymbol{X}}_{2,t})\pi^{j}_{t,1}(\boldsymbol{X}_{1,t})}-\left\\{\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1,t})}-1\right\\}\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1,t})\right.$
$\displaystyle-\left.\left\\{\frac{C^{j}}{\pi^{j}_{t,2}(\overline{\boldsymbol{X}}_{2,t})\pi^{j}_{t,1}(\boldsymbol{X}_{1,t}1)}-\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1,t})}\right\\}\widehat{Q}^{j}_{2,t}(\overline{\boldsymbol{X}}_{2,t})-\theta^{j}\right],$

where $\overline{C}^{j}_{1,t}=I\\{A_{1,t}=d^{j}_{1}(\boldsymbol{X}_{1,t})\\}$,
and for brevity we write
$\widehat{Q}^{j}_{2,t}(\overline{\boldsymbol{X}}_{2,t})=Q_{2}\\{\overline{\boldsymbol{X}}_{2,t},\overline{d}_{2}^{j}(\overline{\boldsymbol{X}}_{2,t};\widehat{\boldsymbol{\beta}}_{2,t}\\}$
and
$\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1,t})=Q^{j}_{1}\\{\boldsymbol{X}_{1,t},d^{j}_{1}(\boldsymbol{X}_{1,t});\widehat{\boldsymbol{\beta}}^{j}_{1,t}\\}$,
which are fitted using $\mathfrak{D}_{t-1}$. To show (i), using the
consistency assumption and rearranging, similar to the argument for the WIPW
estimator, if $W^{A,j}_{t}$ is a function of $\mathfrak{D}_{t-1}$,

$\displaystyle E$
$\displaystyle\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{2,t},\overline{\boldsymbol{A}}_{2,t},Y_{t};\theta^{j})|\mathfrak{D}_{t-1}\\}=W^{A,j}_{t}\left[\vphantom{\frac{C^{j}_{t}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}}E\big{[}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}|\mathfrak{D}_{t-1}\big{]}\right.$
(B.6)
$\displaystyle+E\left(\left.\left[\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}-1\right]\\{Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1})\\}\right|\mathfrak{D}_{t-1}\right)$
(B.7)
$\displaystyle+\left.\left(\left.\left[\frac{C^{j}_{t}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}-\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}\right][Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{2,t}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}]\right|\mathfrak{D}_{t-1}\right)\right].$
(B.8)

Using $\mathcal{W}^{*}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathfrak{D}_{t-1}$, it follows that
$E[\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}|\mathfrak{D}_{t-1}]=E[\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}]=0$,
so that the conditional expectation in (B.6) is zero. Using arguments similar
to those for the WIPW estimator, the conditional expectations in (B.7) and
(B.8) can also be shown to be equal to zero using
$E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})$
and
$E(\overline{C}^{j}_{1,t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\pi_{t,1}^{j}(\boldsymbol{X}_{1})$.
Thus, (i) holds.

To find
$E[\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})\\}^{2}|\mathfrak{D}_{t-1}]$,
first note that the conditional (on $\mathfrak{D}_{t-1}$) expectations of
crossproduct terms in
$\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})\\}^{2}$
are equal to zero by similar arguments. For example, it is straightforward
that

$\displaystyle(W^{A,j}_{t})^{2}E\left(\left.\left[\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}-1\right]\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}\\{Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1})\\}\right|\mathfrak{D}_{t-1}\right)=0$

using $\mathcal{W}^{*}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathfrak{D}_{t-1}$ and
$E(\overline{C}^{j}_{1,t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\pi_{t,1}^{j}(\boldsymbol{X}_{1})$.
Thus,

$\displaystyle E$
$\displaystyle[\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})\\}^{2}|\mathfrak{D}_{t-1}]=(W^{A,j}_{t})^{2}\left[\vphantom{\frac{C^{j}_{t}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}}E[\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}|\mathfrak{D}_{t-1}]\right.$
(B.9)
$\displaystyle+E\left(\left.\left[\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}-1\right]^{2}\\{Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1})\\}^{2}\right|\mathfrak{D}_{t-1}\right)$
(B.10)
$\displaystyle+\left.\left(\left.\left[\frac{C^{j}_{t}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}-\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}\right]^{2}[Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{2,t}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}]^{2}\right|\mathfrak{D}_{t-1}\right)\right].$
(B.11)

In (B.9), using $\mathcal{W}^{*}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathfrak{D}_{t-1}$,
$E[\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}|\mathfrak{D}_{t-1}]=E[\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}]$.
Again using $\mathcal{W}^{*}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\mathfrak{D}_{t-1}$ and the fact that
$E(\overline{C}^{j}_{1,t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\pi_{t,1}^{j}(\boldsymbol{X}_{1})$
implies that

$E\left[\left\\{\left.\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}-1\right\\}^{2}\right|\mathcal{W}^{*},\mathfrak{D}_{t-1}\right]=\frac{1-\pi^{j}_{t,1}(\boldsymbol{X}_{1})}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})},$

so that the conditional expectation in (B.10) is

$E\left[\left.\\{Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1})\\}^{2}\left\\{\frac{1-\pi^{j}_{t,1}(\boldsymbol{X}_{1})}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}\right\\}\right|\mathfrak{D}_{t-1}\right].$

Finally, using
$E(C^{j}_{t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})$
and
$E(\overline{C}^{j}_{1,t}|\mathcal{W}^{*},\mathfrak{D}_{t-1})=\pi_{t,1}^{j}(\boldsymbol{X}_{1})$,
it is straightforward to derive that

$E\left(\left.\left[\frac{C^{j}_{t}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}-\frac{\overline{C}^{j}_{1,t}}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}\right]^{2}\right|\mathcal{W}^{*},\mathfrak{D}_{t-1}\right)=\frac{1-\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi^{j}_{t,1}(\boldsymbol{X}_{1})}.$

Using this result, the conditional expectation in (B.11) is

$E\left(\left.[Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{2,t}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}]^{2}\frac{1-\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi^{j}_{t,1}(\boldsymbol{X}_{1})}\right|\mathfrak{D}_{t-1}\right).$

Substituting the foregoing results in (B.9), (B.10), and (B.11) yields

$\displaystyle E$
$\displaystyle[\\{M_{t}^{j}(\overline{\boldsymbol{X}}_{K,t},\overline{\boldsymbol{A}}_{K,t},Y_{t};\theta^{j})\\}^{2}|\mathfrak{D}_{t-1}]=(W^{A,j}_{t})^{2}\left\\{\vphantom{\frac{C^{j}_{t}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}}E\big{[}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}^{2}\big{]}\right.$
$\displaystyle+E\left[\left.\\{Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1})\\}^{2}\left\\{\frac{1-\pi^{j}_{t,1}(\boldsymbol{X}_{1})}{\pi^{j}_{t,1}(\boldsymbol{X}_{1})}\right\\}\right|\mathfrak{D}_{t-1}\right]$
(B.12)
$\displaystyle+\left.E\left(\left.[Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{2,t}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}]^{2}\left[\frac{1-\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi^{j}_{t,1}(\boldsymbol{X}_{1})}\right]\right|\mathfrak{D}_{t-1}\right)\right\\}.$

Define $\pi^{j}_{t,1}$ and $\pi^{j(s)}_{t,2}$, $s=0,1$, which depend only on
$\mathfrak{D}_{t-1}$ and $j$, and $R^{*}(\mathbf{d}^{j})$ as for the WIPW
estimator, so that
$\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}=I\\{R^{*}(\mathbf{d}^{j})=0\\}\pi^{j(0)}_{t,2}+I\\{R^{*}(\mathbf{d}^{j})=1\\}\pi^{j(1)}_{t,2}$
as in (B.2). Then (B.12) can be rewritten as

$\displaystyle(W^{A,j}_{t})^{2}$
$\displaystyle\left\\{\vphantom{\frac{C^{j}_{t}}{\pi_{t,2}^{j}\\{\overline{\boldsymbol{X}}_{2}^{*}(\mathbf{d}^{j})\\}\pi_{t,1}^{j}(\boldsymbol{X}_{1})}}E\big{[}\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}\big{]}^{2}\right.$
(B.13)
$\displaystyle+\left(\frac{1-\pi^{j}_{t,1}}{\pi^{j}_{t,1}}\right)E\left[\left.\\{Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1})\\}^{2}\right|\mathfrak{D}_{t-1}\right]$
(B.14)
$\displaystyle+\left(\frac{1-\pi^{j(0)}_{t,2}}{\pi^{j(0)}_{t,2}\pi^{j}_{t,1}}\right)E\left(\left.I\\{R^{*}(\mathbf{d}^{j})=0\\}[Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{2,t}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}]^{2}\right|\mathfrak{D}_{t-1}\right)$
(B.15)
$\displaystyle+\left.\left(\frac{1-\pi^{j(1)}_{t,2}}{\pi^{j(1)}_{t,2}\pi^{j}_{t,1}}\right)E\left(\left.I\\{R^{*}(\mathbf{d}^{j})=1\\}[Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{2,t}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}]^{2}\right|\mathfrak{D}_{t-1}\right)\right\\}.$
(B.16)

Write the expectations in (B.13)–(B.16) respectively as
$\nu^{j}=E[\\{Y^{*}(\mathbf{d}^{j})-\theta^{j}\\}]^{2}$,
$\nu^{j}_{1}=E[\\{Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{1,t}(\boldsymbol{X}_{1})\\}^{2}|\mathfrak{D}_{t-1}]$,
and

$\nu^{j(s)}_{2}=E\left(\left.I\\{R^{*}(\mathbf{d}^{j})=1\\}[Y^{*}(\mathbf{d}^{j})-\widehat{Q}^{j}_{2,t}\\{\overline{\boldsymbol{X}}_{k}^{*}(\mathbf{d}^{j})\\}]^{2}\right|\mathfrak{D}_{t-1}\right),\,\,\,s=0,1,$

so that the variance (B.12) can be written succinctly as

$\displaystyle(W^{A,j}_{t})^{2}\left\\{\nu^{j}+\nu^{j}_{1}\left(\frac{1-\pi^{j}_{t,1}}{\pi^{j}_{t,1}}\right)+\nu^{j(0)}_{2}\left(\frac{1-\pi^{j(0)}_{t,2}}{\pi^{j(0)}_{t,2}\pi^{j}_{t,1}}\right)+\nu^{j(1)}_{2}\left(\frac{1-\pi^{j(1)}_{t,2}}{\pi^{j(1)}_{t,2}\pi^{j}_{t,1}}\right)\right\\}.$
(B.17)

Thus, as for the WIPW estimator, we wish to choose $W^{A,j}_{t}$ so that
(B.17) is a constant depending on $j$ but not $\mathfrak{D}_{t-1}$. Denoting
this constant as $\Xi^{A,j}$, the weights should be chosen so that

$W^{A,j}_{t}=(\Xi^{A,j})^{1/2}\left\\{\nu_{j}+\nu^{j}_{j}\left(\frac{1-\pi^{j}_{t,1}}{\pi^{j}_{t,1}}\right)+\nu^{j(0)}_{2}\left(\frac{1-\pi^{j(0)}_{t,2}}{\pi^{j(0)}_{t,2}\pi^{j}_{t,1}}\right)+\nu^{j(1)}_{2}\left(\frac{1-\pi^{j(1)}_{t,2}}{\pi^{j(1)}_{t,2}\pi^{j}_{t,1}}\right)\right\\}^{-1/2}.$
(B.18)

With a burn-in period of length $t^{*}$ weeks of nonadaptive randomization
with $W^{A,j}_{t}\equiv 1$, from (B.18), it must be that

$\Xi^{A,j}=\nu^{j}+\nu^{j}_{1}\left(\frac{1-\pi^{j}_{t^{*},1}}{\pi^{j}_{t^{*},1}}\right)+\nu^{j(0)}_{2}\left(\frac{1-\pi^{j(0)}_{t^{*},2}}{\pi^{j(0)}_{t^{*},2}\pi^{j}_{t^{*},1}}\right)+\nu^{j(1)}_{2}\left(\frac{1-\pi^{j(1)}_{t^{*},2}}{\pi^{j(1)}_{t,2}\pi^{j}_{t^{*},1}}\right),$
(B.19)

where as before the randomization probabilities during the burn-in period are
fixed, so that $\pi^{j(s)}_{t^{*},2}$, $s=0,1,$ and $\pi_{t^{*},1}^{j}$ are
known constants not depending on the data. Thus, if $\nu^{j}$, $\nu^{j}_{1}$,
and $\nu^{j(s)}_{2}$, $s=0,1$, were known, so that $\Xi^{A,j}$ is known, the
weights for $t\geq t^{*}+1$ should be chosen as in (B.18) with $\Xi^{A,j}$ as
in (B.19).

Analogous to the approach for the WIPW estimator, the proposed implementation
scheme is as follows. Obtain weights by estimating $\Xi^{A,j}$ by an estimator
$\widehat{\Xi}^{A,j}_{t^{*}}$ based on the data from a burn-in period of
length $t^{*}$ weeks with nonadaptive randomization with $W^{A,j}_{t}\equiv
1$, $t\leq t^{*}$, and then going forward treat the estimate as fixed and
known. Letting $\widehat{\nu}^{j}_{t}$, $\widehat{\nu}^{j}_{t,1}$, and
$\widehat{\nu}^{j(s)}_{t,2}$, $s=0,1$, be estimators for $\nu^{j}$,
$\nu^{j}_{1}$, and $\nu^{j(s)}_{2}$, $s=0,1$, at any time $t$ based on
$\mathfrak{D}_{t-1}$, define

$\widehat{\Xi}^{A,j}_{t}=\widehat{\nu}^{j}_{t}+\widehat{\nu}^{j}_{t,1}\left(\frac{1-\pi^{j}_{t,1}}{\pi^{j}_{t,1}}\right)+\widehat{\nu}^{j(0)}_{t,2}\left(\frac{1-\pi^{j(0)}_{t,2}}{\pi^{j(0)}_{t^{*},2}\pi^{j}_{t,1}}\right)+\widehat{\nu}^{j(1)}_{t,2}\left(\frac{1-\pi^{j(1)}_{t,2}}{\pi^{j(1)}_{t,2}\pi^{j}_{t,1}}\right).$
(B.20)

Then (B.19) implies that the estimator $\widehat{\Xi}^{A,j}_{t^{*}}$ for
$\Xi^{A,j}$ based on the burn-in data is

$\widehat{\Xi}^{A,j}_{t^{*}}=\widehat{\nu}^{j}_{t^{*}}+\nu^{j}_{1}\left(\frac{1-\pi^{j}_{t^{*},1}}{\pi^{j}_{t^{*},1}}\right)+\widehat{\nu}^{j(0)}_{t^{*},2}\left(\frac{1-\pi^{j(0)}_{t^{*},2}}{\pi^{j(0)}_{t^{*},2}\pi^{j}_{t^{*},1}}\right)+\widehat{\nu}^{j(1)}_{t^{*},2}\left(\frac{1-\pi^{j(1)}_{t^{*},2}}{\pi^{j(1)}_{t,2}\pi^{j}_{t^{*},1}}\right).$

Then as for the WIPW estimator, treating $\widehat{\Xi}^{A,j}_{t^{*}}$ as
fixed and known, for $t\geq t^{*}+1$, from (B.18), take
$W^{A,j}_{t}=(\widehat{\Xi}^{A,j}_{t^{*}}/\widehat{\Xi}^{A,j}_{t})^{1/2}$.
Because $\widehat{\Xi}^{A,j}_{t^{*}}$ depends only on