# Unveiling Bias in Sequential Decision Making:
A Causal Inference Approach for Stochastic Service Systems

Juan Camilo David
Department of Industrial and Systems Engineering
University of Wisconsin-Madison
Madison, WI
<EMAIL_ADDRESS>
& Amy Cochran
Department of Population Health Sciences
and Department of Mathematics
University of Wisconsin-Madison
Madison, WI
<EMAIL_ADDRESS>
& Gabriel Zayas-Cabán
Department of Industrial and Systems Engineering
University of Wisconsin-Madison
Madison, WI
<EMAIL_ADDRESS>

###### Abstract

In many stochastic service systems, decision-makers find themselves making a
sequence of decisions, with the number of decisions being unpredictable. To
enhance these decisions, it is crucial to uncover the causal impact these
decisions have through careful analysis of observational data from the system.
However, these decisions are not made independently, as they are shaped by
previous decisions and outcomes. This phenomenon is called _sequential bias_
and violates a key assumption in causal inference that one person’s decision
does not interfere with the potential outcomes of another. To address this
issue, we establish a connection between sequential bias and the subfield of
causal inference known as dynamic treatment regimes. We expand these
frameworks to account for the random number of decisions by modeling the
decision-making process as a marked point process. Consequently, we can define
and identify causal effects to quantify sequential bias. Moreover, we propose
estimators and explore their properties, including double robustness and
semiparametric efficiency. In a case study of 27,831 encounters with a large
academic emergency department, we use our approach to demonstrate that the
decision to route a patient to an area for low acuity patients has a
significant impact on the care of future patients.

_Keywords_ Sequential bias, causal inference, marked point process, semi-
parametric efficiency theory, doubly robust estimation.

##  1 Introduction

Sequential decision-making problems are pervasive and thus the subject of
extensive analysis in Operations Research and Management Sciences, Computer
Science, Statistics, and other disciplines. These problems are primarily
motivated by the necessity of enhancing a sequence of decisions for a _single_
unit, such as a person or single job. One example is medical decision-making,
which searches for strategies that optimize outcomes for an individual
patient. This is often achieved using Markov decision processes (Alagoz et al.
2010, Steimle and Denton 2017). Another example is dynamic treatment regimes
(DTRs), which aims to estimate the causal impact of sequences of treatment
decisions compared to a baseline strategy. These estimates are then used to
find a sequence of decisions that enhances outcomes for each patient. This
approach has considerable attention in the statistical literature (Robins
1997, Murphy 2003, Chakraborty and Moodie 2013, Tsiatis et al. 2019).

Instead of examining multiple decisions for a single unit, we explore a
scenario where a single decision-maker navigates a sequence of decisions over
time for _multiple_ units. Our objective is not to identify optimal strategies
for such scenarios, but rather to assess and estimate the direct or causal
influence that one unit’s decision has on future unit’s decisions and
outcomes. Inspired by the psychology literature, we refer to this influence as
_sequential bias_ (Yu and Cohen 2008).

Service systems have provided many compelling examples of sequential bias. The
authors’ inspiration derives from their research of split-flow models in
emergency departments (EDs), where the triage nurse is replaced by a physician
(Zayas-Cabán et al. 2016, David et al. 2022). As part of split flow, the
physician at triage decides between routing patients to a vertical area for
low-acuity patients or a traditional room. In resource-constrained settings
like the ED, physicians must weigh the care of the current patient against the
needs of future patients. This prompts the question of whether a physician’s
routing decision for one patient influences the care of future patients.

Another example, with which researchers studying service systems are likely
familiar, is the organ transplantation process. Led by the Organ Procurement
and Transplantation Network, the transplantation process involves a list of
patients awaiting organs, prioritized based on availability, type, location,
and health status. When an organ becomes available, potential recipients are
identified according to compatibility and urgency. A patient’s transplant
team, acting as decision-makers, assess organ compatibility and faces a time-
sensitive binary choice of accepting or rejecting the organ. Due to the
limited supply of organs, this sequential decision-making process is
susceptible to sequential bias, wherein the current decision carries life-and-
death consequences for future patients.

Sequential bias is relevant in other domains. In many selection processes, for
example, evaluators assess candidates in a sequential manner to determine
their suitability for a role, position, or opportunity. These processes can
include grading students, judging sporting competitions, determining criminal
sentencing, or approving loans, among others. The prior evaluations of
candidates can impact subsequent scores assigned by evaluators, potentially
introducing bias. The order in which candidates are evaluated plays a crucial
role as it can affect the outcomes and favor certain individuals while
disadvantaging others (Chen et al. 2016, Goldbach et al. 2022).

A common thread runs through all these examples, allowing for a single
conceptual framework. A decision-maker encounters various jobs within a
stochastic service system and must make decisions for each job sequentially
over time. The number of jobs being intervened upon is random and may be
beyond the immediate control of the decision-maker. Each job shares a common
set of binary intervention options, and the decision for each job bears
consequences for the job itself. These decisions are determined based on
contextual information and the history of prior interventions, contexts, and
outcomes. Additionally, the decision-maker lacks information on future
contexts or outcomes when making the current decision. The objective is to
quantify sequential bias using observational data from the system. This
entails defining, identifying, and estimating the causal effect that a
decision for the current job has on the decisions and outcomes of future jobs.
By accurately quantifying this bias, we unlock opportunities to improve
decision-making in stochastic service systems through stochastic modeling and
optimization.

Performing causal inference within this context presents two issues. First,
there is the issue of interference, where the decision for one job influences
decisions and outcomes of future jobs (Cox 1958). Interference violates a core
assumption in causal inference known as the Stable Unit Treatment Value
Assumption (SUTVA) (Rubin 1980). An immediate implication is that the
potential decisions and outcomes of a future job is dependent on the entire
history of past interventions. The collection of these histories grows
exponentially in the number of jobs, making it impractical to contrast all
possible sequences of interventions. The DTR literature also encounters this
exponential growth, offering us ideas for focusing on specific causal
contrasts (Guo et al. 2021). Causal blip and excursion effects are examples of
such focused approaches, as they provide more manageable linear growth in the
number of decisions (Robins 1997, Boruvka et al. 2018, Qian et al. 2021).

Second, there is the issue of a random number of variables, further
complicated by the possible dependence of this random number on the history of
prior interventions, contexts, and outcomes. This brings into question on how
causal frameworks, such as Pearl’s do-calculus (Pearl 2000), Neyman-Rubin
potential outcomes (Neyman 1923, Rubin 1974), or specific frameworks for DTRs
(Murphy 2003), can be readily applied. These frameworks are typically intended
for scenarios where there are a fixed and finite number of random variables.
One option would be to set an upper bound on the random number of jobs.
However, it is often preferred to avoid such a restrictive assumption in
stochastic service systems. Consequently, the random number of jobs leads us
to, in effect, deal with an infinite number of random variables, presenting
challenges across various aspects. These challenges include ensuring the
existence of probability distributions, handling integrability, and carefully
managing the interchangeability of summation and differentiation in this
infinite context.

In this paper, we contribute a solution to these issues as follows. We define
a causal model to capture the sequential decision-making scenario and the
random number of jobs. To address the varying number of variables, we
artificially extend the variables to represent a marked point process
(Jacobsen and Gani 2006). This extension results in an infinite number of
random variables. Subsequently, we utilize Markov kernels to construct a
probability distribution for the marked point process and to represent
hypothetical scenarios were decisions to be modified. This construction aligns
with the causal frameworks proposed by Richardson and Robins (2013) and
Malinsky et al. (2019), which unified Pearl’s do-calculus (Pearl 2000) with
the Neyman-Rubin potential outcomes framework. It also facilitates the use of
the single-world intervention graph (SWIG) and encompasses common causal
inference assumptions such as consistency and sequential ignorability. For
clarity, we illustrate the causal model with directed acyclic graphs
restricted to a finite number of variables.

The subsequent step is defining and identifying causal effects. These causal
effects are defined in a manner that closely resembles the causal contrasts
utilized by Boruvka et al. (2018) and Qian et al. (2021), while also
accounting for the adjustments needed to accommodate the random number of
original variables. We introduce two types of causal effects: the lag effect
and the marginalized lag effect, where the latter is a dimensionally reduced
version of the former. A crucial consideration for these effects is the need
to condition on post-intervention variables. This is to ensure that we
exclusively measure the causal impact on future jobs that actually exist,
given that the number of jobs is random. For insightful discussion on
conditioning on post-intervention variables, refer to Pearl (2015).
Importantly, we demonstrate that our causal effects can be nonparametrically
identified.

The next step is estimation, for which we propose an estimating equation
approach with careful consideration for integrability. We prove our estimator
is doubly-robust and apply standard statistical arguments to establish its
consistency and asymptotic normality (Van der Vaart 2000). Furthermore, using
techniques from semiparametric theory (Tsiatis 2006), we derive the efficient
score for a specific case of our lag effect, which supports our choice of
estimator. Last, we apply our estimator to electronic healthcare records
(n=27,831) from a large academic hospital that operates a split-flow model in
the ED. We investigate the causal impact that routing the current patient to a
vertical area over a traditional room has on the subsequent patient.

The remainder of this article is structured as follows. Section 2 summarizes
relevant studies from the literature, encompassing topics such as DTRs and
sequential randomized trials, interference in causal inference, and causal
reasoning for stochastic processes. It also highlights the contribution of
this study in relation to the surveyed literature. Section 3 introduces the
causal model employed in this paper, utilizing marked point processes, and
defines the causal effects of interest. Additionally, it presents conditions
under which these causal effects can be identified from observational data.
Moving forward, Section 4 presents an estimator for these causal effects,
accompanied by derivations of its asymptotic properties such as double
robustness and its efficient score. Section 5 applies the estimation procedure
to the split-flow example discussed earlier, aiming to determine the impact of
prior routing decisions on the subsequent patient’s outcomes. Section 6
concludes our study.

##  2 Literature Review

###  2.1 Dynamic Treatment Regimes and Sequential Randomized Trials

The current study relates to the subfield of causal inference that explores
sequential decision-making for individual units. Specifically, it closely
aligns with DTRs, sequential multiple assignment randomized trials (SMARTs),
and their high-frequency counterparts known as micro-randomized trials.

A DTR is a specified sequence of treatment decision rules that determines how
to adapt the delivery of treatments over time in response to an individual’s
changing health and other time-varying contextual factors. A DTR typically
consists of multiple stages, where each stage considers an individual’s
medical history and current health information to recommend the next
treatment. Extensive research has focused on determining the optimal DTR
(Murphy 2003, Laber et al. 2014, Li et al. 2023) and estimating the causal
effects of specific treatment regimes relative to an alternative treatment
regime (Robins 1986, 1997, Murphy 2003, Chakraborty and Murphy 2014). One of
the key challenges addressed in this literature is the exponential growth in
the total number of treatment regimes as the number of treatment decisions
increases. For instance, binary decisions at four stages results in a total of
$2^{4}$ possible treatment sequences. As a consequence, the amount of
available data for comparing two treatment sequences significantly decreases
as more decision stages are added. This reduction in data introduces greater
uncertainty in the comparisons and optimization of DTRs.

When faced with an overwhelming number of treatment sequences, a commonly
employed approach is to focus on DTRs that can be represented using observed
treatment decisions. These DTRs are typically characterized by a linear
increase in their complexity with respect to treatment stages. These concepts
are captured in modeling frameworks like the structural nested mean model
(Robins 1994) and quantified through measures such as blip effects (Robins
1997, Wang and Yin 2020) and causal excursion effects Qian et al. (2021), Shi
et al. (2022), among others (Guo et al. 2021). To illustrate, a blip effect
can compare treatment regimes that initially align with the observed treatment
until a certain stage, after which each regime discontinues the treatment. The
sole distinction between these regimes lies in the timing of when the
treatment is discontinued. In general, these causal effects are typically
described as the average response of an individual at a specific stage,
considering two different treatment regimes, while conditioning on their
previous treatment history and time-dependent factors. The DTR literature
employs standard assumptions for causal inference, such as positivity,
consistency, and conditional exchangeability. These assumptions play a crucial
role in non-parametrically identifying the causal effects and enable us to
perform causal inference even in non-experimental settings, such as when
working with observational data.

Related to DTRs are experimental designs known as SMARTs and also referred to
as alternative designed randomized trials by Robins (1986). These designs
involve repeatedly assigning participants to different treatment conditions at
multiple stages, considering their individual characteristics and the
historical information of covariates and treatments (Lavori and Dawson 2000,
2004). These trials help inform DTRs, providing at evidence of a randomized
control trial with the advantage of tailoring treatment decisions more
efficiently based on a patient’s current and past context. For example, if a
patient does not respond to a treatment, they can be re-randomized to
alternative treatments in order to find a more effective approach. The seminal
work of Murphy (2005) presents a robust framework for designing SMARTs.
Moreover, the statistical methods of a SMART can be adapted to non-
experimental contexts with a few key distinctions. Specifically, assumptions
regarding certain conditional independence relations are necessary, and the
randomization probabilities of treatment assignments typically need to be
estimated from data. To safeguard inferences against incorrect assumptions,
robust methods are recommended (Guo et al. 2021).

Similar to SMARTs, micro-randomized trials involve repeatedly randomizing
participants to different treatment conditions over time (Klasnja et al.
2015). These trials build on the causal inference concepts for analyzing DTRs
from longitudinal data, introduced in Robins (1997) and Gill and Robins
(2001). In a micro-randomized trial, randomization occurs frequently, often
multiple times a day. This high-frequency randomization allows researchers to
adapt interventions to an individual’s rapidly changing context, such as
variations in stress or mood. Micro-randomized trials are especially suitable
for evaluating mobile health interventions (Klasnja et al. 2015, Boruvka et
al. 2018, Qian et al. 2021) and developing just-in-time adaptive interventions
(Nahum-Shani et al. 2017). As an example, an author of this study conducted
micro-randomized trials to evaluate a mobile version of a therapy called
acceptance and commitment therapy (Cochran et al. 2023, Thomas et al. 2023).
Extensive efforts have been made to devise statistical methods for micro-
randomized trials (Klasnja et al. 2015, Boruvka et al. 2018, Qian et al.
2021).

When designing DTRs, SMARTs, and micro-randomized trials, decision epochs
represent predetermined points in time when treatment decisions are made for
each individual. Notably, decision epochs have a finite and fixed number in
these designs. For example, a SMART frequently uses two epochs: a trial’s
onset and a later time when treatment response is evaluated. Meanwhile, micro-
randomized trials incorporate multiple decision epochs throughout each day
over a specified duration, such as fixed morning and evening time windows over
six weeks in the aforementioned trials (Cochran et al. 2023, Thomas et al.
2023). Consequently, existing statistical methods have historically overlooked
the added complexity of a random number of decision epochs.

By contrast, stochastic service systems encounter a random number of decision
epochs when performing sequential decision-making. This variability can arise
from the unpredictable number of arrivals to the service system and the random
duration required to service jobs within the system. In a stochastic service
system, the number of decision epochs is not only random but can also be
influenced by the system’s history. Decisions may be restricted to specific
time periods (e.g., 9am to 5pm) or restricted to a fixed number of jobs
receiving a higher level of resources. Due to the system’s stochastic nature,
these constraints result in a random number of decisions.

###  2.2 Interference in Causal Inference

A stochastic service system involves a sequence of decisions for different
jobs or individuals, which differs from DTRs, SMARTs, and micro-randomized
trials where sequential decisions are often centered around the same
individual. This presents a challenge because most literature on causal
inference relies on SUTVA (Rubin 1980, 2005). Part of SUTVA is the no
interference assumption (Cox 1958), which says that the potential outcomes of
a specific unit are unaffected by the decisions on other units. The sequential
decision processes analyzed in this work violate the no interference
assumption, as the decisions made for previous jobs can impact the potential
outcomes of future jobs. Consequently, we must appropriately account for
interference when identifying and estimating causal effects in the context of
sequential decision-making within a stochastic service system.

Early approaches to interference involved dividing units into equal-sized and
non-overlapping blocks, allowing interference within blocks but not across
them, known as partial interference (Hong and Raudenbush 2006, Sobel 2006,
Hudgens and Halloran 2008). Recent advancements have allowed for more
arbitrary patterns of interference (Verbitsky-Savitz and Raudenbush 2012,
Sofrygin and van der Laan 2017). These approaches permit interference among
subsets of units based on spatial proximity, social ties, and other measures
of proximity. An example of modeling general interference without partial
interference is presented by Aronow and Samii (2017), where a function is
introduced to model exposure using aggregated data on the number of exposed
neighbors. This approach is flexible as it allows an arbitrary number of
neighbors. However, it has the limitation of not allowing other types of
interactions, such as contagion, defined as the outcome of a unit causing the
outcome of a different unit. More recently, Tchetgen Tchetgen et al. (2021)
proposed a fairly flexible method for inferring causal inferences from complex
networks, and Zhang et al. (2022) build a model where specific forms of
interactions (e.g., spillovers, contagion) can be explicitly incorporated. The
study by Zhang et al. (2022) also develops an algorithm to quantify
interaction bias and conditions under which the bias may be ignored. However,
the framework is restricted to causal linear models.

Our view is that interference within a wide range of stochastic service
systems is notably less complex compared to the arbitrary interference
highlighted earlier. There are two reasons for this. First, stochastic service
systems have a sequential nature, meaning that future events cannot causally
influence past events. This sequentiality leads to a progressive accumulation
of information over time, akin to a filtration in the context of a stochastic
process. Interestingly, this sequential interference gives rise to similar
relationships of conditional independence between observed variables as seen
in DTRs, SMARTs, and micro-randomized trials. As a result, the approaches used
to handle conditional independence in the latter can be applied to address
interference in stochastic service systems.

Second, we can anticipate partial interference in numerous real-world service
systems, which can be attributed to two underlying mechanisms. One mechanism
is the operation of service systems within fixed time windows, such as 9am to
5pm. This holds for the ED split-flow model that served as the inspiration for
our research. It is reasonable to assume that interference occurs among jobs
within the same time window but not between jobs across different time
windows. A second mechanism is the use of distinct decision-makers in a
service system. There are situations where it is justifiable to assume that
interference occurs among jobs handled by the same decision-maker but not
between jobs assigned to different decision-makers. Thus, the concept of
partial interference, as discussed in previous studies (Hong and Raudenbush
2006, Sobel 2006, Hudgens and Halloran 2008), can be applied. Specifically, we
adopt an assumption of independent data panels, whereby interference occurs
between jobs within each panel, but not between jobs from different panels.

###  2.3 Stochastic Processes and Causal Inference

Existing literature on causal inference for stochastic systems remains
relatively scarce, but is increasingly gaining recognition (Didelez 2008,
Røysland 2011, 2012, Didelez 2015, Gao et al. 2021). Significantly, many
foundational elements commonly employed in causal inference in other contexts,
such as causal models and graphical representations, face significant
challenges when applied to stochastic systems. Early representative work is
the study by Didelez (2008) which proposes a new class of graphical models
capturing the dependence structure of events that occur in time via marked
point processes. Although our research does not directly employ these
graphical models, their work inspired our adoption of marked point processes.
Another distinction is our interest in identification and estimation of causal
effects using observational data, as opposed to graphical analyses.

##  3 Causal Framework

###  3.1 Notation

We use the following notation:

  * •

Capital letters for random variables (e.g., $Z$).

  * •

Greek letters for parameters (e.g., $\theta,\xi,\eta,\alpha,\beta$).

  * •

${\cal X}^{(Z)}$ for the target space of a random variable $Z$.

  * •

Lower case letters for observations of a random variable (e.g., $z\in{\cal
X}^{(Z)}$).

  * •

Subscripts for vectors or sequences, e.g., $Z_{1:k}$ is $(Z_{1},\ldots,Z_{k})$
and $Z_{B}$ is $(Z_{i_{1}},\ldots,Z_{i_{k}})$ when
$B=\\{i_{1},\ldots,i_{k}\\}\subseteq\mathop{\rm I\kern-1.99997ptN}\nolimits$.

  * •

Apostrophe for transpose, e.g. $Z^{\prime}$ is the transpose of $Z$.

  * •

$\|Z\|$ for the $L_{2}$-norm of $Z$.

  * •

$\nabla_{\alpha}$ for partial differentiation with respect to $\alpha$.

  * •

$\mathop{\rm I\kern-1.99997ptP}\nolimits^{(Z)}$ for a Markov kernel from some
source ${\cal X}^{(W)}$ to a target ${\cal X}^{(Z)}$, which are regular
conditional distributions that generalize the transition matrix for Markov
chains to general state spaces.

We adopt the convention that any such Markov kernel is a probability
distribution on ${\cal X}^{(Z)}$ when the source ${\cal X}^{(W)}$ is
$\emptyset$ and can be viewed as a Markov kernel with source ${\cal
X}^{(W,V)}$ that is insensitive to some of its arguments, viz,

$\mathop{\rm I\kern-1.99997ptP}\nolimits^{(Z)}(\cdot|w,v)=\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(Z)}(\cdot|w),\quad\forall\,(w,v)\in{\cal
X}^{(W,V)}.$

In addition, we adopt the convention that $Z_{1:k}$ and ${\cal X}^{Z_{1:k}}$
are empty when $k=0$ and that ${\cal X}^{(Z_{1:2})}$ represents the product
space ${\cal X}^{(Z_{1})}\times{\cal X}^{(Z_{2})}$ associated with $Z_{1:2}$
(similarly for sequences).

###  3.2 Setting

In our stochastic service system, jobs (such as customers, patients, or units)
are assumed to arrive one at a time and are processed individually. Each job
is associated with specific characteristics like age, sex, race, chief
complaint, comorbidities, and the current capacity of the system. A decision-
maker makes a binary decision for each job, which can be influenced by the
job’s characteristics and by information on jobs in the past. The decision
leads to an outcome for the job, such as the length of time they stay in the
system or the number of tests performed. In this system, the decision-maker
processes a _random_ number of jobs in any given shift or day, which is called
a panel. It is assumed that multiple panels are associated with the service
system.

We collect data on a panel using several variables. Let $K\geq 1$ be the
random number of jobs in a panel, and let $k$ denote the index of a job within
a panel ($k=1,\ldots,K$). Each job $k$ is associated with three random
variables:

  * •

$X_{k}\in\mathbb{R}^{l}$ (a vector of characteristics associated with the
job),

  * •

$A_{k}\in\\{0,1\\}$ (a binary intervention), and

  * •

$Y_{k}\in\mathbb{R}$ (an outcome of interest).

We assume that the set of variables from a single panel are mutually
independent across panels. The history of all three variables up to and
including job $k$ is denoted by $H_{k}=(X_{1:k},A_{1:k},Y_{1:k})$. We assume
an ordering of the jobs, which means that a decision $A_{k}$ can only depend
on information available immediately prior to the decision, including the
job’s present characteristics $X_{k}$, and the history $H_{k-1}$ consisting of
past job characteristics $X_{1:k-1}$, decisions $A_{1:k-1}$, and outcomes
$Y_{1:k-1}$. This assumption is expressed more formally later.

###  3.3 A Marked Point Process

One of the foundations of causal inference is a causal model describing a
probability model and causal relationships between variables. The variables
are usually fixed and finite and either discrete or absolutely continuous.
This poses a challenge for this paper’s interest in a random number of jobs.
To overcome this challenge, we propose using marked point processes (MPPs) as
a probability model for the variables. This requires extending the definition
of the random variables.

The first step in this extension is to attach a time $T_{k}$ to each job $k$,
constrained so that

$0<T_{1}<T_{2}\ldots<T_{K}<\infty,$

where $T_{k}$ could be the observed or unobserved time at which the system
first processes job $k$. To define times for all $k\in N$, we let
$T_{k}=\infty$ whenever $k>K$. The $T_{k}$ can now comprise the _times_ of a
marked point process (MPP), which must be an increasing (infinite) sequence
taking values in $(0,\infty]$. The original random number $K$ of jobs can be
recovered from the times $T_{k}$ by taking the supremum of the set
$\\{k:T_{k}<\infty\\}$.

The _marks_ of a marked point process also need to be defined for all
$k\in\mathop{\rm I\kern-1.99997ptN}\nolimits$. This is done by taking
$(X_{k},A_{k},Y_{k})$ to be the marks of the process when $k\leq K$, and
introducing an irrelevant mark $\Delta$ to extend the definition of these
variables to all $k\in N$. Specifically, $X_{k}$, $A_{k}$, and $Y_{k}$ can
take irrelevant marks $\Delta$ to signify that the corresponding time $T_{k}$
is infinite, or equivalently, that the panel size $K$ is less than $k$. The
irrelevant mark is useful in subsequent derivations to remind us that it only
makes sense to manipulate the marks numerically when $k\leq K$. Moreover, the
original variables can be recovered perfectly from the $T_{k}$ and
$(X_{k},A_{k},Y_{k})$.

The second step of this extension is to define appropriate measurable spaces
for the times and marks. For times up to a finite value $n$, we use the set

${\cal X}^{(T_{1:n})}=\big{\\{}t_{1:n}\in\mathbb{R}^{n}:0<t_{1}\leq\ldots\leq
t_{n};\,t_{k}<t_{k+1}\text{ if }t_{k}<\infty\big{\\}},$

which contains strictly increasing times while finite, and increasing times in
general. Similarly, for an infinite sequence of times, we define the analogous
set ${\cal X}^{(T_{1:\infty})}$. For the marks, we use the sets

$\displaystyle{\cal X}^{(X_{k})}=\mathbb{R}^{l}\cup{\Delta},\quad{\cal
X}^{(A_{k})}=\\{0,1\\}\cup{\Delta};\quad{\cal
X}^{(Y_{k})}=\mathbb{R}\cup{\Delta},$

which consist of marks that are either real-valued or the irrelevant mark
$\Delta$. By taking products of these sets, additional sets such as ${\cal
X}^{(X_{1:k})}$, ${\cal X}^{(X_{1:\infty})}$, ${\cal X}^{(H_{k})}$, and ${\cal
X}^{(H_{\infty})}$ can be defined. These spaces are equipped with appropriate
Borel $\sigma$-algebras and product Borel $\sigma$-algebras, as described in
Jacobsen and Gani (2006).

By extending the variable definition, we can model the distribution
$\mathop{\rm I\kern-1.99997ptP}\nolimits$ for
$(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$. We
can induce $\mathop{\rm I\kern-1.99997ptP}\nolimits$ by constructing random
variables iteratively from Markov kernels, which follows from the Ionescu
Tulcea theorem (c.f. Theorem B.3.5. in Hinderer et al. (2016) and Theorem A.1
in Appendix). Section 3.4 provides a more detailed specification of the Markov
kernels. Once these kernels are properly defined, we can view the collection
$(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$ as a
_marked point process_. This is a stochastic process that resides in ${\cal
X}^{(T_{1:\infty})}\times{\cal X}^{(H_{\infty})},$ for which the $T_{k}$ tend
to infinity, and $X_{k}$, $A_{k}$, and $Y_{k}$ are the irrelevant mark if and
only if $T_{k}$ is infinite.

###  3.4 Causal Model

To place our random variables within a causal inference framework, we follow
the approach of Richardson and Robins (2013) and Malinsky et al. (2019), with
one modification: we replace the use of conditional densities with Markov
kernels. This allows us to handle an infinite collection of general random
variables. A causal model is defined by specifying two objects for each random
variable $Z$. The first object is a set of variables, denoted by
$\mathop{{\rm{Pa}}}\nolimits(Z)$, of _parents_ or _direct causes_ of $Z$. In
general, the parents are those variables that, when intervened upon, directly
affect $Z$. We make the following assumptions about the parents:

###### Assumption 1

(Causal model - parents) For each $k\in\mathop{\rm
I\kern-1.99997ptN}\nolimits$, let

  * •

$\mathop{{\rm{Pa}}}\nolimits(T_{k})=(T_{1:k-1},H_{k-1})$

  * •

$\mathop{{\rm{Pa}}}\nolimits(X_{k})=(T_{1:k},H_{k-1})$

  * •

$\mathop{{\rm{Pa}}}\nolimits(A_{k})=(X_{k},H_{k-1})$

  * •

$\mathop{{\rm{Pa}}}\nolimits(Y_{k})=(X_{k},A_{k},H_{k-1}).$

These sets give rise to a directed acyclic graph (DAG), where the nodes
represent variables and a directed edge exists from one variable to another if
the former is a parent of the latter.

To allow for greater generality, we tried not to impose any restrictions on
which past observations could be direct causes of future observations. For
instance, the parents of $T_{k}$ include all past observations: $T_{1:k-1}$
and $H_{k-1}$, and the same can be said for $X_{k}$. However, in order to
identify the causal effects of interest, we had to impose some restrictions.
For example, the parents of $Y_{k}$ and $A_{k}$ do not include the possibly
unobserved times $T_{1:k}$.

The DAG can be used to visualize the sets of parents. Although we cannot
visualize the entire DAG for our variables, we can visualize the DAG
restricted to a finite set. Figure 1 shows the DAG restricted to
$(T_{k},X_{k},A_{k},Y_{k})_{k=1,2}$. It is worth noting that every DAG can be
topologically ordered, which means that the nodes can be ordered in a way that
ensures no node comes before its parents.

Figure 1: DAG induced by assumed causal model restricted to
$(T_{k},X_{k},A_{k},Y_{k})_{k=1,2}$ when $K\geq 2$.
$T_{1}$$X_{1}$$A_{1}$$Y_{1}$$T_{2}$$X_{2}$$A_{2}$$Y_{2}$

The second object is a Markov kernel from ${\cal
X}^{(\mathop{{\rm{Pa}}}\nolimits(Z))}$ to ${\cal X}^{(Z)}$ for each $Z$. These
Markov kernels tell us how to construct a variable $Z$ in
$(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$ from
its parents. They also tell how to construct new variables were we to modify
the value of its parents. To ensure we can easily move between each set of
variables, it is desirable to define the same probability measure for both.
Thus, in addition to assuming the existence of Markov kernels, it benefits us
to regard the original and new variables as arising from the same collection
of independent variables, often referred to as _exogenous_ variables. These
considerations motivate defining Markov kernels as arising from a probability
distribution for an exogenous variable and a deterministic assignment:

###### Assumption 2

(Causal model - Markov kernels) Assuming we have parents given by Assumption
1, each variable $Z$ in the collection
$(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$ is
associated with a Markov kernel $\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(Z)}$ from ${\cal
X}^{(\mathop{{\rm{Pa}}}\nolimits(Z))}$ to ${\cal X}^{(Z)}$. We assume this
Markov kernel can be expressed as:

$\displaystyle\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(Z)}(B\,|\,\mathop{{\rm{Pa}}}\nolimits_{Z})=\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(\varepsilon_{Z})}(\\{w\in\Omega^{(Z)}\,|\,f^{(Z)}(\mathop{{\rm{Pa}}}\nolimits_{Z},w)\in
B\\}),$

for some probability distribution $\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(\varepsilon_{Z})}$ on $(\Omega^{(Z)},{\cal
E}^{(Z)})$ and a measurable function $f^{(Z)}$ from ${\cal
X}^{(\mathop{{\rm{Pa}}}\nolimits(Z))}\times\Omega^{(Z)}$ to ${\cal X}^{(Z)}.$

To model $(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}$ as an MPP, we require further conditions on the
Markov kernels ensuring that, almost surely, the times are strictly increasing
when finite and increasing otherwise, and the marks take the irrelevant mark
exactly when the corresponding time is infinite. These assumptions are stated
precisely in the Appendix (Assumption A.1). Our causal model, now comprised of
sets of parents and Markov kernels, gives us a probability model for the
variables we are interested in:

###### Assumption 3

(Causal model - distribution) Suppose that we have parents and kernels as
given by Assumptions 1,2, and A.1. We obtain the distribution $\mathop{\rm
I\kern-1.99997ptP}\nolimits$ of a marked point process for
$(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$ in
two steps:

  * •

Obtain a distribution $\mathop{\rm I\kern-1.99997ptP}\nolimits$ on
$(\varepsilon_{T_{k}},\varepsilon_{X_{k}},\varepsilon_{A_{k}},\varepsilon_{Y_{k}})_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}$ using the Ionescu Tulcea theorem with respect to
their respective sequence of probability distributions.

  * •

Assign iteratively
$Z=f^{(\varepsilon_{Z})}(\mathop{{\rm{Pa}}}\nolimits(Z),\varepsilon_{Z})$ for
each variable in ${(T_{k},X_{k},A_{k},Y_{k})}_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}$

Finally, assume $K=\sup\\{k:T_{k}<\infty\\}$ under $\mathop{\rm
I\kern-1.99997ptP}\nolimits$ is finite almost surely.

It is important to mention that the ordering in which variables are assigned
reflects a topological ordering of our variables that is compatible with our
DAG. Specifically, each variable $Z$ must be defined after its parents are. In
addition, we could have applied the Ionescu Tulcea theorem directly to the
Markov kernels $\mathop{\rm I\kern-1.99997ptP}\nolimits^{(Z)}$ to arrive at
the same distribution $\mathop{\rm I\kern-1.99997ptP}\nolimits$. However, as
mentioned previously, it would be advantageous to fix the exogenous variables
while manipulating the values of parents, so all our variables are defined
with the same distribution $\mathop{\rm I\kern-1.99997ptP}\nolimits$.

###  3.5 Potential Outcomes

Markov kernels describe the conditional distribution of a variable as a
function of its parents. Crucially, by replacing the parents in the Markov
kernel with fixed values, we can examine what would have happened under
different scenarios. Specifically, we focus on examining what would have
occurred if we had fixed the values of $A_{k}$:

######  Definition 1 (Potential outcomes)

Suppose we have $\mathop{\rm I\kern-1.99997ptP}\nolimits$ on
${(\varepsilon_{T_{k}},\varepsilon_{X_{k}},\varepsilon_{A_{k}},\varepsilon_{Y_{k}})}_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}$ from Assumption 3. We define the potential
outcomes for $B\subseteq\mathop{\rm I\kern-1.99997ptN}\nolimits$ and
$a\in{\cal X}^{(A_{1:\infty})}$,

$(T_{k}(a_{B}),X_{k}(a_{B}),A_{k}(a_{B}),Y_{k}(a_{B}))_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits},$

iteratively according to
$Z(a_{B})=f^{(\varepsilon_{Z})}(\mathop{{\rm{Pa}}}\nolimits[a_{B}](Z),\varepsilon_{Z}),$
where we obtain $\mathop{{\rm{Pa}}}\nolimits[a_{B}](Z)$ by replacing each
variable $W$ in ${\mathop{{\rm{Pa}}}\nolimits}(Z)$ either with $a_{k}$ if $W$
is $A_{k}$ and $k\in B$ or with $W(a_{B})$ otherwise. The sets
$\mathop{{\rm{Pa}}}\nolimits[a_{B}](Z)$ give rise to a single world
intervention graph (SWIG), where the nodes represent potential outcomes and
the $a_{k}$, $k\in B$, and a directed edge exists from one variable to
$Z(a_{B})$ if the former is in the set
$\mathop{{\rm{Pa}}}\nolimits[a_{B}](Z)$.

To put it another way, if the variable $A_{k}$ for $k\in B$ was originally
used to construct a subsequent variable, the potential outcome definition
dictates we pass forward $a_{k}$ instead of $A_{k}(a_{B})$. Even though
$A_{k}(a_{B})$ does not get passed forward, we still construct $A_{k}(a_{B}),$
which distinguishes the causal inference framework of Richardson and Robins
(2013) from others. Moreover, while other frameworks may make certain
assumptions about potential outcomes, assumptions like consistency end up
being a natural consequence (c.f., Malinsky et al. (2019)):

###### Proposition 1

(Consistency) Suppose we have potential outcomes defined for $B\cup
C\subseteq\mathop{\rm I\kern-1.99997ptN}\nolimits$ and $C\subseteq\mathop{\rm
I\kern-1.99997ptN}\nolimits$ with $B\cap C=\emptyset$ and $a\in{\cal
X}^{(A_{1:\infty})}$. Among events for which $A_{B}(a_{B\cup C})=a_{B}$, the
potential outcomes $Z(a_{B})$ and $Z(a_{B\cup C})$ are equal.

To clarify, consistency implies that potential outcomes $Z(a_{B\cup C})$
matches $Z(a_{C})$ when the potential outcome $A_{B}(a_{C})$ matches $a_{B}$.
In essence, we need consistency so that $Y_{k}(a_{B})$ has the same
distribution as the original variable $Y_{k}$ when the intervention $A_{B}$
matches the value $a_{B}$ that we have fixed for the intervention. This
enables us to substitute $Y_{k}(a_{B})$ for $Y_{k}$ when we condition on
$A_{B}=a_{B}$. The proof of consistency follows by noting that, under the
conditions of the proposition, $Z(a_{B})$ and $Z(a_{B\cup C})$ are expressed
identically in terms of the exogenous variables $\varepsilon^{(W)}$ and the
assignment functions $f^{(W)}$ for $W$ in
${(T_{k},X_{k},A_{k},Y_{k})}_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}.$
This identity allows us to claim the potential outcomes are equal _everywhere_
, as opposed to the weaker _almost surely_ or _in distribution_.

The definition of potential outcomes leads us to another important concept
called causal irrelevance, which tells us when we can drop the potential
outcome notation. The SWIG induced by the potential outcome definition is the
most convenient way to express causal irrelevance:

###### Proposition 2

(Causal irrelevance) Suppose we have potential outcomes defined for $B\cup
C\subseteq\mathop{\rm I\kern-1.99997ptN}\nolimits$ and $C\subseteq\mathop{\rm
I\kern-1.99997ptN}\nolimits$ with $B\cap C=\emptyset$ and $a\in{\cal
X}^{(A_{1:\infty})}$. If $Z(a_{B\cup C})$ is not a descendent of any $a_{k}$
with $k\in C$ in our SWIG, then the potential outcomes $Z(a_{B})$ and
$Z(a_{B\cup C})$ are equal.

In other words, if a potential outcome is not a descendant in the SWIG of an
intervened value $a_{k}$, for $k\in B$, then the potential outcome is said to
be causally irrelevant to the choice of $a_{k}$. In such cases, the value of
the potential outcome does not change, whether or not we had passed $a_{k}$
forward at all. For example, since $T_{1}$ has no parents, then $T_{1}(a_{B})$
is always equal to $T_{1}$. The proof of causal irrelevance is similar to that
of consistency, i.e. under the conditions of the proposition, $Z(a_{B})$ and
$Z(a_{B\cup C})$ are expressed identically in terms of the exogenous variables
and the assignment functions.

After applying causal irrelevance, we visualize the SWIG restricted to
$(T_{k}(a_{1}),X_{k}(a_{1}),A_{k}(a_{1}),Y_{k}(a_{1}))_{k=1,2}$ with
$B=\\{1\\}$ in Figure 2. To obtain the SWIG from the original DAG, we split
each node corresponding to $A_{k}$ for $k\in B$ into two nodes, one inheriting
all incoming edges of $A_{k}$ and the other inheriting all outgoing edges.

Figure 2: SWIG restricted to
$(T_{k}(a_{1}),X_{k}(a_{1}),A_{k}(a_{1}),Y_{k}(a_{1}))_{k=1,2}$ when $K\geq
2$.$T_{1}$$X_{1}$$A_{1}$ $a_{1}$
$Y_{1}(a_{1})$$T_{2}(a_{1})$$X_{2}(a_{1})$$A_{2}(a_{1})$$Y_{2}(a_{1})$

One of the key advantages of using SWIGs is that they uncover conditional
independence relationships among potential outcomes. These relationships are a
result of the SWIG being _Markov compatible_ with the potential outcome
distribution, which then implies local and global Markov properties. The local
Markov property, for example, says that a variable is independent of its non-
descendants when conditioning on its parents. The global Markov property
characterizes conditional independence relationships arising from a property
of the Markov compatible graph, known as _d-separation_. While we will not
delve into these concepts further in this paper, we will use a specific
application of the global Markov property of the SWIG:

###### Proposition 3

(Sequential ignorability) Suppose we have potential outcomes defined for
$B=\\{k\\}$ and $a\in{\cal X}^{(A_{1:\infty})}$. Then, $Y_{k+\ell}(a_{k})$ is
independent of $A_{k}$ conditional on $T_{k}(a_{k})<\infty$ and any vector
$R_{k}(a_{k})$ of variables that includes $H_{k-1}$, $X_{k}$, and a (possibly
empty) subset of variables in
$\mathop{{\rm{Pa}}}\nolimits[a_{k}](Y_{k+\ell})$.

To demonstrate the above is an implementation of the global Markov property,
we consider all paths in the SWIG that connect $A_{k}$ to $Y_{k+\ell}(a_{k})$.
Because $A_{k}$ has no outgoing edges (the SWIG shifts these edges to
$a_{k}$), then each of these paths must go through the parents of $A_{k}$.
These parents are exactly $X_{k}$ and $H_{k-1}$, which are included in
$R_{k}$. In addition, it is important to note that the initial parent
encountered by every path cannot be a collider. This is because $A_{k}$ is not
considered a parent to any of its own parents. These conditions imply that
$R_{k}(a_{k})$ d-separates $A_{k}$ from $Y_{k+\ell}(a_{k})$ in the SWIG, which
in turn gives Proposition 3. We refer the reader to Richardson and Robins
(2013) and Malinsky et al. (2019) for more on Markov properties of SWIGs.

###  3.6 Lag Effect

Now that we have our causal model in place and our potential outcomes defined,
we can focus on what we want to learn from the data. Our aim is to understand
the potential outcome of a future job if the current job is subjected to a
specific intervention. Let $\ell>0$ denote the time lag between the current
job being intervened upon and a future job that might be affected by it. We
refer to job $k+\ell$ as the future job and job $k$ as the current job. We can
consider the potential outcomes as defined in Definition 1 with $B=\\{k\\}$
and $a_{k}\in\\{0,1\\}$ to capture this idea. According to causal irrelevance
(Proposition 2), we have $T_{i}(a_{k})=T_{i}$, $X_{i}(a_{k})=X_{i}$ and
$A_{i}(a_{k})=A_{i}$ when $i\leq k$; $Y_{i}(a_{k})=Y_{i}$ when $i<k$.
Otherwise, we cannot drop the potential outcome notation. Provided it is not
the irrelevant mark, we are interested in $Y_{k+\ell}(a_{k})$, which
represents the outcome for the future job we would have observed if the
current job was forced to take intervention $a_{k}$. Although we cannot learn
about individual potential outcomes $Y_{k+\ell}(a_{k})$, we can estimate their
expectations in various situations:

######  Definition 2 (Lag effects)

Suppose we have potential outcomes given by Definition 1 for $B=\\{k\\}$ and
$a_{k}=1$ and separately for $a_{k}=0$. Let $R_{k}$ be any vector of variables
that includes $H_{k-1}$, $X_{k}$, and any (possibly empty) subset of variables
in $\mathop{{\rm{Pa}}}\nolimits(Y_{k+\ell})\setminus A_{k}$. Introduce for
$k,\ell\in\mathop{\rm I\kern-1.99997ptN}\nolimits$,

$\displaystyle\zeta(k,\ell,r)=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}(1)|R_{k}(1)=r,T_{k+\ell}(1)<\infty\right]-$
$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}(0)|R_{k}(0)=r,T_{k+\ell}(0)<\infty\right].$

to represent the average causal effect on job $k+\ell$ after intervening on
job $k$ among specific events for which $T_{k+\ell}(a_{k})<\infty$ and
$R_{k}(a_{k})=r$ under the respective intervention.

In the above definition, we compute expectation with respect to the
probability distribution $\mathop{\rm I\kern-1.99997ptP}\nolimits$, while
holding the interventions of the first $k-1$ jobs at their original values and
varying only the intervention for job $k$. The lag effect may depend on the
index $k$ of the current job and the index of $k+\ell$ of the future job. To
avoid computing the irrelevant mark, we condition on
$T_{k+\ell}(a_{k})<\infty$.

It is important to understand the causal query being asked with the lag
effect, as it is subtle. We expose job $k$ to intervention $a_{k}=1$ and
average the resulting value of $Y_{k+\ell}$ over all events which
_consequently_ attain a certain level for $R_{k}$ (possibly post-intervention)
and yield observations at time period $k+\ell$. This scenario is then
contrasted to exposing job $k$ to intervention $a_{k}=0$ instead and averaging
the resulting value of $Y_{k+\ell}$ over all events which consequently attain
the same level for $R_{k}$ and yield observations at time period $k+\ell$. Our
causal query is also distinctly different than what it would be were we to
use:

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}(1)|R_{k}=r,T_{k+\ell}<\infty\right]-$
$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}(0)|R_{k}=r,T_{k+\ell}<\infty\right].$

Above asks what would happen under different interventions among events who
_actually_ attain a certain level of $R_{k}$ and yield observations at time
period $k+\ell$ under the realized intervention. While neither are
conventional, Pearl has a nice discussion on distinguishing between these two
causal queries (Pearl 2015) and why someone might condition on post-
intervention variables.

For us, we have chosen this lag effect based on practical considerations. One
important consideration is that we need to condition on at least one post-
intervention event, specifically the event where the induced $Y_{k+\ell}$ is
observed. This is necessary to avoid computing expectations of irrelevant
marks, and is unavoidable when the number of jobs, $K$, is random and
influenced by interventions. This situation is common in many stochastic
service systems, such as those that process jobs at different rates depending
on prior interventions and in a fixed time frame (e.g., 9a to 5p).

Another factor to consider is that although the latter causal query, which
conditions on realized values, may be more straightforward to comprehend, it
is generally not identifiable. The latter query is equivalent to conditioning
on colliders (e.g., $T_{k+\ell}$) between the intervention $A_{k}$ and the
potential outcome $Y_{k+\ell}(a_{k})$, thereby opening up a path to transmit
non-causal associations from $A_{k}$ to $Y_{k+\ell}(a_{k})$. The only way to
close the path is to condition on potential outcomes related to post-
intervention variables (such as $T_{k+\ell}(a_{k})$). The two causal queries
would coincide if the conditioning event
$\\{R_{k}(a_{k})=r,T_{k+\ell}(a_{k})<\infty\\}$ remains invariant to $a_{k}$.
For example, if $R_{k}$ only includes pre-intervention variables such as
$H_{k-1}$ and $X_{k}$, and the intervention has no impact on whether the
subsequent $Y_{k+\ell}$ is observed, then
$\\{R_{k}(a_{k})=r,T_{k+\ell}(a_{k})<\infty\\}$ could be invariant to $a_{k}$.
Yet, as noted above, we anticipate that the latter condition may be too
restrictive for many stochastic systems.

A final consideration was what variables should we allow in the conditioning
set. We choose to require $H_{k-1}$ and $X_{k}$ in $R_{k}$ for reasons that we
will discuss in the next section, as they are relevant to identifying the
causal effect. We allow for variables that occur post-intervention, but prior
to $Y_{k+\ell}$, so we later can better model the variation in $Y_{k+\ell}$.
For example, the future job’s characteristics $X_{k+\ell}$ may have a larger
influence on the outcome $Y_{k+\ell}$ than the current job’s characteristics
$X_{k}$. We also consider these post-intervention variables to focus on the
direct effects of the intervention, which fix the induced value of these
variables to a certain level. For instance, we may wish to determine if
directing a patient to vertical flow in the ED would speed up service for a
future job, even were the induced congestion levels held constant.

Because the lag effect is high-dimensional, we may find it too difficult to
specify the right functional form to express the influence of $r$ on the lag
effect $\zeta(k,\ell,r)$. In these cases, it may be beneficial to use
marginalization to reduce the dimension of the lag effect and then express the
functional form for this lower-dimensional effect. One way to accomplish this
is with the following:

######  Definition 3 (Marginalized lag effects)

Consider the lag effect $\zeta(k,\ell,r)$ as it is given in Definition 2. Let
$S_{k}$ be any (possibly empty) subset of $R_{k}$. Introduce

$\displaystyle\zeta_{\rm{marg}}(k,\ell,s)=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\zeta(k,\ell,R_{k})|S_{k}=s,T_{k+\ell}<\infty\right].$

to represent the lag effect marginalized over $R_{k}$ conditional on $S_{k}=s$
and $T_{k+\ell}<\infty$.

If $S_{k}=R_{k}$, then the lag effect and the marginalized lag effect are
equivalent. Our choice of the marginalized lag effect is also guided by
practical considerations. In this case, our concern is not about
identification. Identifying the marginalized lag effect can be achieved once
the lag effect is identified. However, it is possible to condition on
potential outcomes instead of observed variables. This would result in an
identifiable effect, which can be calculated through a parametric g-formula
(Robins 1986). However, this approach requires a full specification of the
Markov kernels, which we find too restrictive. As we saw before, if the
conditioning event $\\{R_{k}(a_{k})=r,T_{k+\ell}(a_{k})<\infty\\}$ remains
invariant to $a_{k}$, then the two strategies coincide, and the marginalized
lag effect is simply:

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}(1)|S_{k}=s,T_{k+\ell}<\infty\right]-$
$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}(0)|S_{k}=s,T_{k+\ell}<\infty\right].$

###  3.7 Identification

The primary challenge in causal inference is that we can only observe one of
the potential outcomes, namely either $Y_{k+\ell}(1)$ or $Y_{k+\ell}(0)$. As a
result, we are generally unable to estimate the lag effect (as defined in
Definition 2). We need to extrapolate from our observations of $Y_{k+\ell}(1)$
when the intervention is assigned to $A_{k}=1$ in order to make inferences
about $Y_{k+\ell}(1)$ when the intervention is assigned to $A_{k}=0$. For such
extrapolation to be valid, we require the assumption of positivity in addition
to consistency, causal irrelevance, and sequential ignorability (as outlined
in Propositions 1–3). This assumption states that every job has a chance of
being assigned to the intervention:

###### Assumption 4

(Positivity) Consider $R_{k}$ as defined for the lag effect (Definition 2).
Assume

$\mathop{\rm
I\kern-1.99997ptE}\nolimits[A_{k}|R_{k},T_{k+\ell}<\infty]\in(0,1)$

almost surely.

With positivity and the other properties, we are able to identify the lag
effect, which is to say that we can re-express the lag effect equivalently in
terms of observed variables:

###### Lemma 1

(Non-parametric identification) Consider the lag effect $\zeta(k,\ell,r)$ as
it is given in Definition 2. Assuming positivity (Assumption 4), then the lag
effect $\zeta(k,\ell,r)$ is equivalent to

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}\,|\,A_{k}=1,R_{k}=r,T_{k+\ell}<\infty]-\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}\,|\,A_{k}=0,R_{k}=r,T_{k+\ell}<\infty].$

#### Proof.

Sequential ignorability (Proposition 3) and positivity (Assumption 4) implies
that

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}(a_{k})\,|\,R_{k}(a_{k})=r,T_{k+\ell}(a_{k})<\infty]=\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}(a_{k})\,|\,A_{k}=a_{k},R_{k}(a_{k})=r,T_{k+\ell}(a_{k})<\infty].$

Consistency (Proposition 1) simplifies the last expression to

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}\,|\,A_{k}=a_{k},R_{k}=r,T_{k+\ell}<\infty].$

Applying these expressions to $a_{k}=1$ and $a_{k}=0$, we arrive at

$\displaystyle\zeta(k,\ell,r)$ $\displaystyle=\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}(1)\,|\,R_{k}(1)=r,T_{k+\ell}(1)<\infty]-\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}(0)\,|\,R_{k}(0)=r,T_{k+\ell}(0)<\infty]$
$\displaystyle=\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}\,|\,A_{k}=1,R_{k}=r,T_{k+\ell}<\infty]-\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}\,|\,A_{k}=0,R_{k}=r,T_{k+\ell}<\infty],$

completing the proof.

##  4 Estimation

###  4.1 Estimating Equation

Our goal is to estimate marginalized lag effects $\zeta_{\rm{marg}}(k,\ell,s)$
using estimating equations for a fixed $\ell\in\mathop{\rm
I\kern-1.99997ptN}\nolimits$. Suppose we have collected data on $n$ panels.
This data is assumed to be comprised of $n$ independent and identically
distributed (iid) realizations of our marked point process
$(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$
modeled with distribution $\mathop{\rm I\kern-1.99997ptP}\nolimits$. We do not
assume $T_{k}$ is observed or use irrelevant marks for estimation, so that
observations are restricted to the set
$(X_{k},A_{k},Y_{k})_{k\in\\{1,\ldots,K\\}}$. Let $\mathop{\rm
I\kern-1.99997ptE}\nolimits$ denote expectation with respect to $\mathop{\rm
I\kern-1.99997ptP}\nolimits$ and $\mathop{\rm I\kern-1.99997ptE}\nolimits_{n}$
denote the empirical average with respect to $n$ iid realizations of our
marked point process.

The primary estimating function we work with is of the form

$\displaystyle\mathop{\rm I\kern-1.99997ptE}\nolimits_{n}[U]=0,$ (1)

where

$\displaystyle U=\sum_{k=1}^{K-\ell}U_{k},\qquad\text{and}\qquad
U_{k}=W_{k}\left(Y_{k+\ell}-g(R_{k})^{\prime}\alpha-
A_{k}f(S_{k})^{\prime}\beta\right)\begin{bmatrix}g(R_{k})\\\
A_{k}f(S_{k})\end{bmatrix}$

for parameters $\alpha,\beta,\xi,\eta$ that we wish to estimate and suitable
functions $g$ and $f$ that depend on $\xi$ and $\eta$. We search for
parameters $\theta:=(\xi,\eta,\alpha,\beta)$ that solve Equation (1) in some
compact subset $\Theta$ of $\mathbb{R}^{m}$. To simplify notation, we dropped
the dependence of functions on parameters. The variable $U_{k}$ can be
decomposed into several parts:

  * •

$W_{k}$ (a weight for mitigating bias and emphasizing certain data points)

  * •

$Y_{k+l}$ (outcome of future job)

  * •

$g(R_{k})^{\prime}\alpha$ (a linear working model for the baseline conditional
mean of $Y_{k+\ell}$)

  * •

$f(S_{k})^{\prime}\beta$ (a linear working model for the marginalized lag
effect)

  * •

$\begin{bmatrix}g(R_{k})\\\ A_{k}f(S_{k})\end{bmatrix}$ (derivative of the
working mean model $g(R_{k})^{\prime}\alpha+A_{k}f(S_{k})^{\prime}\beta$ with
respect to $\alpha$, $\beta$)

In particular, the parameter $\beta$ is the inferential target, as it
contributes to the marginalized linear working model for the lag effect. If
$K$ were deterministic, then the estimating equation (Equation 1) has several
precedents. It is akin to a generalized estimating equation (GEE), in which
data are weighted, clustered into panels, and given an independent working
correlation matrix. It is also akin to the g-estimation approach of
Vansteelandt and Sjolander (2016) for structural nested mean models and the
weighted and centered approach of Boruvka et al. (2018).

While not the inferential target, the estimating equation depends on
parameters $\xi$ and $\eta$. These parameters specify a working model for
conditional probabilities of job $k$’s intervention assignment:

$\displaystyle q_{k}(S_{k};\xi)$ $\displaystyle\approx\mathop{\rm
I\kern-1.99997ptE}\nolimits[A_{k}\,|\,S_{k},T_{k+\ell}<\infty]$ $\displaystyle
p_{k}(R_{k};\eta)$ $\displaystyle\approx\mathop{\rm
I\kern-1.99997ptE}\nolimits[A_{k}\,|\,R_{k},T_{k+\ell}<\infty]$

We assume $\xi$ is also estimated using an estimating equation of the form

$\displaystyle\mathop{\rm I\kern-1.99997ptE}\nolimits_{n}[L]=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}\left[\sum_{k=1}^{K-\ell}L_{k}\right]=0,$

for variables $L_{k}\mathbbm{1}_{\\{k\leq K\\}}$ that depends only on $\xi$
and the observations of $(X_{k},A_{k},Y_{k})_{k\in\\{1,\ldots,K\\}}$.
Similarly, we assume $\eta$ is also estimated using an estimating equation of
the form

$\displaystyle\mathop{\rm I\kern-1.99997ptE}\nolimits_{n}[M]=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}\left[\sum_{k=1}^{K-\ell}M_{k}\right]=0$

for variables $M_{k}\mathbbm{1}_{\\{k\leq K\\}}$ that depends only on $\eta$
and the observations of $(X_{k},A_{k},Y_{k})_{k\in\\{1,\ldots,K\\}}$. Most
likely, $\xi$ and $\eta$ would be estimated using logistic or log-linear
regression. For example, if we used logistic regression and assumed a linear
model in $S_{k}$ on the log-odds scale with $\xi$ as the regression
coefficients, our estimating equation might use:

$L=\sum_{k=1}^{K-\ell}(A_{k}-\mathrm{logit}^{-1}([1,S_{k}]^{\prime}\xi))\begin{bmatrix}1\\\
S_{k}\end{bmatrix}.$

By stacking the estimating equations together,

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}\left\\{\begin{bmatrix}U\\\ L\\\
M\end{bmatrix}\right\\}=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}\left\\{\sum_{k=1}^{K-\ell}\begin{bmatrix}U_{k}\\\
L_{k}\\\ M_{k}\end{bmatrix}\right\\}=0,$ (2)

estimation can be viewed as searching for $\alpha$, $\beta$, $\xi$, and $\eta$
that solve the stacked estimating equation. If we define $U_{\rm{stacked}}$ as
the concatenated vector of $U$, $L$, and $M$, then in order for Equation (2)
to be meaningful, we are assuming that $U_{\rm{stacked}}$ is integrable. Since
these variables are sums of a random number of variables, we impose mild
conditions on $K$ and $U_{k}$, $L_{k}$, and $M_{k}$ to recover integrability
(see Assumption A.2 in Appendix) to recover integrability:

######  Lemma 2 (Integrability)

Under the conditions on $U_{k}$, $L_{k}$, $M_{k}$ and $K$ in Assumption A.2,
Equation 2 is well-defined in the sense that, for each $\theta\in\Theta$,
$U_{\rm{stacked}}$ is integrable and has finite expectation.

The assumptions can be found in Appendix A, and the proof in Appendix B. The
broad arc is to replace the random sum in the definition of $U_{\rm{stacked}}$
with an infinite sum and then exchange the expectation and the infinite sum.
The conditions in Assumption A.2 are then used to ensure that the final sum
can be bounded above by the expectation of $K$ up to a constant.

###  4.2 Estimating Procedure

The first step of estimation would be to use the observations of
$(X_{k},A_{k},Y_{k})_{k\in\\{1,\ldots,K\\}}$ to recover an estimate $\xi_{n}$
and $\eta_{n}$ that solves the estimating equations $\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}[U]=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}[V]=0$. For example, we might use logistic
regression to model the above conditional probabilities on the log odds scale
as a linear function of the $S_{k}$ or $R_{k}$. In such a case, the parameters
$\xi_{n}$ and $\eta_{n}$ would be the estimated regression coefficients. In
special circumstances, the intervention $A_{k}$ was randomized according to a
known probability. Consequently, $q_{k}(S_{k};\xi)$ and $p_{k}(R_{k};\eta)$
would be exact.

Once $\xi$ and $\eta$ are estimated, we turn to specifying $W_{k}$, $f$, and
$g$. The weight term $W_{k}$ is taken to be the ratio of the two conditional
probability models:

$\displaystyle
W_{k}=A_{k}\frac{q_{k}(S_{k};\xi_{n})}{p_{k}(R_{k};\eta_{n})}+(1-A_{k})\frac{1-q_{k}(S_{k};\xi_{n})}{1-p_{k}(R_{k};\eta_{n})}.$
(3)

Naturally, we are assuming that we are not dividing by zero, otherwise the
weights would be ill-defined. This assumption is the empirical equivalent of
our positivity assumption (Assumption 4). These weights resemble (stabilized)
inverse probability weights, in which a job is weighted inversely according to
the probability of receiving their realized intervention assignment (Horvitz
and Thompson 1952). Jobs that are more likely to receive their realized
intervention assignment are subsequently down-weighted relative to those that
are less likely. Given that certain jobs may be predisposed to be assigned a
particular intervention, the weights create a pseudo-population of jobs in
each intervention group that better reflects all jobs, not just those in the
given intervention group. In the special case when $S_{k}=R_{k}$ or when
$p_{k}(R_{k};\eta_{n})=q_{k}(S_{k};\xi_{n})$, then the weights $W_{k}$ are
just $1$.

Meanwhile, we specify $f(S_{k})$ to be any vector-valued function of $S_{k}$.
It reflects our assumptions about how we think the marginalized lag effect
might vary as a function of $S_{k}$, i.e.

$\displaystyle f(S_{k})^{\prime}\beta\approx\xi_{\rm{marg}}(k,\ell,S_{k})$

For example, we might think the marginalized lag effect of assigning a patient
to vertical flow is quadratic in a patient’s age. If $S_{k}$ were the current
patient’s age, then we might chose:

$\displaystyle
f(S_{k})=\begin{bmatrix}1&S_{k}&S_{k}^{2}\end{bmatrix}^{\prime}$

We can similarly specify the baseline term $g(R_{k})$ to be any vector-valued
function of $R_{k}$ with one caveat. It must contain the variable
$q_{k}(S_{k};\xi)f(S_{k})$. This caveat is so that the working mean model:

$\displaystyle g(R_{k})^{\prime}\alpha+A_{k}f(S_{k})^{\prime}\beta,$

upon a re-definition of $\alpha$, includes the model in which the $A_{k}$ is
centered:

$\displaystyle
g(R_{k})^{\prime}\alpha+(A_{k}-q_{k}(S_{k};\xi_{n}))f(S_{k})^{\prime}\beta$

Including $q_{k}(S_{k};\xi)f(S_{k})$ in $g(R_{k})$ follows the strategy of
Vansteelandt and Sjolander (2016) of making the mean model more general while
still, practically, centering the intervention assignments. Centering is
recommended for reasons we discuss later and is recommended in Boruvka et al.
(2018).

The last step of estimation is to search for a solution $\alpha_{n}$ and
$\beta_{n}$ to the estimating equation (Equation 1). Given the similarity of
Equation (1) to GEE, these parameters (though not necessarily their standard
errors) can be estimated using standard GEE software, provided the working
correlation matrix is specified to be independent. We next study the
properties of the resulting estimator $\beta_{n}$.

###  4.3 Consistency and Asymptotic Normality

Our estimator is a type of estimator called a Z-estimator (Z for “zero"), as
it involves searching for roots of an estimating equation. Asymptotic
properties of Z-estimators have been extensively studied, so that we can
invoke standard arguments. In an effort to be self-contained, we sketch the
arguments in Van der Vaart (2000). We start with the following definition:

######  Definition 4 (Estimator)

Define an estimator $\theta_{n}=(\xi_{n},\eta_{n},\alpha_{n},\beta_{n})$ to be
a solution in a set $\Theta\subseteq\mathbb{R}^{m}$ to the (stacked)
estimating equation (Equation 2), should such a solution exist. Define
$\theta_{\infty}=(\xi_{\infty},\eta_{\infty},\alpha_{\infty},\beta_{\infty})$
to be a solution in $\Theta$ to the asymptotic version of the estimating
equation:

$\mathop{\rm I\kern-1.99997ptE}\nolimits\left[U_{\rm{stacked}}\right]=0,$

should such a solution exist.

The first property of our estimators is consistency:

######  Theorem 1 (Consistency)

Under Assumption A.3, the estimator $\theta_{n}$ in Definition 4 is consistent
with respect to $\theta_{\infty}$:

$\theta_{n}\stackrel{{\scriptstyle P}}{{\rightarrow}}\theta_{\infty}.$

#### Proof.

The proof follows from Theorem 5.7 in Van der Vaart (2000) if two conditions
are met. Take $M_{n}(\theta)=-\begin{Vmatrix}\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}[U_{\rm{stacked}}]\end{Vmatrix}$ and
$M(\theta)=-\begin{Vmatrix}\mathop{\rm
I\kern-1.99997ptE}\nolimits[U_{\rm{stacked}}]\end{Vmatrix}$. The first
condition is the uniform convergence of $M_{n}(\theta)$ to $M(\theta)$ in
probability over $\Theta$. Given that $\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}[U_{\rm{stacked}}]$ is the empirical average
over iid samples of $U_{\rm{stacked}}$, this uniform convergence is an
application of a uniform law of large numbers and follows from several
conditions, including compactness of $\Theta$ in $\mathbb{R}^{m}$, finiteness
of $\mathop{\rm I\kern-1.99997ptE}\nolimits[U_{\rm{stacked}}]$ (Lemma 2), and
the mapping $\theta\rightarrow U_{\rm{stacked}}$ being almost surely
continuous and dominated by an appropriate function (see Lemma 1 in Tauchen
(1985)). The second condition is called well-separatedness and requires that
the supremum of $M(\theta)$ over all $\Theta$ outside of any neighborhood
around a unique maximizer $\theta_{\infty}$ is positive. Well-separatedness
follows immediately from the uniqueness of the maximizer and the almost sure
continuity of the mapping $\theta\rightarrow U_{\rm{stacked}}$ over compact
$\Theta$ (Assumptions A.3). The proof of Theorem 5.7 in Van der Vaart (2000)
uses the uniform convergence condition to argue that
$M(\theta_{n})\stackrel{{\scriptstyle P}}{{\rightarrow}}M(\theta_{\infty})=0$.
Well-separatedness then ensures that the only way $M(\theta_{n})$ converges in
probability to $0$ would be if $\theta_{n}$ were converging in probability to
$\theta$.

The assumptions required for consistency (Assumption A.3) in common choices of
$L$ and $M$ (such as logistic or log-linear models) are relatively mild.
Continuity of $\theta\rightarrow U_{\rm{stacked}}$ is reasonable, since $U$ is
linear in $\alpha$ and $\beta$ and smooth in $\xi$, and $L$ and $M$ is
constant in $\alpha$ and $\beta$ and smooth in $\xi$ and $\eta$. The main
concern is that $U$ involves division by $p_{k}(R_{k};\eta)$ or
$1-p_{k}(R_{k};\eta)$, which could prevent continuity over $\Theta$ if
$p_{k}(R_{k};\eta)$ is not bounded away from $0$ and $1$. For a function to
dominate the mapping $\theta\rightarrow U_{\rm{stacked}}$, $K$ and variables
contributing to $U_{k}$ cannot be too large. This is comparable to the
conditions for integrability (Assumption A.2). In practice, these variables
are usually bounded. Similarly, parameters are usually bounded in practice, so
compactness of $\Theta$ is usually not an issue.

The second property of our estimators is $\sqrt{n}$-consistency and asymptotic
normality:

######  Theorem 2 (Asymptotic normality)

Under Assumption A.4, the estimator $\theta_{n}$ in Definition 4 is
$\sqrt{n}$-consistent with respect to $\theta_{\infty}$:

$\sqrt{n}(\theta_{n}-\theta_{\infty})=-\mathbb{B}^{-1}\sqrt{n}E_{n}\left[U_{\rm{stacked}}\big{|}_{\theta=\theta_{\infty}}\right]+o_{p}(1)$

and asymptotically normal with mean $\theta_{\infty}$ and variance
$\mathbb{B}^{-1}\mathbb{C}(\mathbb{B}^{-1})^{\prime}$, where

$\displaystyle\mathbb{B}=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\nabla_{\theta}U_{\rm{stacked}}\big{|}_{\theta=\theta_{\infty}}\right];\qquad\mathbb{C}=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[U_{\rm{stacked}}U_{\rm{stacked}}^{\prime}\big{|}_{\theta=\theta_{\infty}}\right].$

#### Proof.

The proof follows from several results in Van der Vaart (2000), including
Theorem 5.21, Lemma 19.24, and Example 19.7. Let $\Psi(\theta)=\mathop{\rm
I\kern-1.99997ptE}\nolimits[U_{\rm{stacked}}]$ and
$\Psi_{n}(\theta)=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}[U_{\rm{stacked}}]$, such that
$\Psi(\theta_{\infty})=\Psi_{n}(\theta_{n})=0$. We can obtain an approximation
for $\sqrt{n}\Psi(\theta_{n})$ from $\Psi(\theta_{\infty})=0$ and the
differentiability of $\Psi(\theta)$ at $\theta_{\infty}$:

$\displaystyle\sqrt{n}\Psi(\theta_{n})=\sqrt{n}\,\mathbb{B}\cdot(\theta_{n}-\theta_{\infty})+\sqrt{n}\begin{Vmatrix}\theta_{n}-\theta_{\infty}\end{Vmatrix}o_{P}\left(1\right).$

Meanwhile, the central limit theorem allows us to conclude that
$\sqrt{n}\Psi_{n}(\theta_{\infty})$ is asymptotically normal with mean zero
and covariance $\mathbb{C}.$ The challenge is to connect the two
aforementioned results through empirical process theory, which necessitates
demonstrating that
$\sqrt{n}\Psi_{n}(\theta_{\infty})+\sqrt{n}\Psi(\theta_{n})\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$. This connection relies on the compactness of $\Theta$
and the (almost sure) Lipschitz continuity of $\theta\rightarrow
U_{\rm{stacked}}$, which provides uniform convergence in distribution over
$\Theta$ of the empirical process
$\sqrt{n}\left(\Psi_{n}(\theta)-\Psi(\theta)\right)$. Thus, we can deduce:

$\displaystyle\sqrt{n}\,\Psi(\theta_{n})=\sqrt{n}\,\mathbb{B}\cdot(\theta_{n}-\theta_{\infty})+\sqrt{n}\,\begin{Vmatrix}\theta_{n}-\theta_{\infty}\end{Vmatrix}\,o_{P}\left(1\right)=-\sqrt{n}\,\Psi_{n}(\theta_{\infty})+o_{P}(1).$

Last, we argue
$\sqrt{n}\,\begin{Vmatrix}\theta_{n}-\theta_{\infty}\end{Vmatrix}\,o_{P}\left(1\right)=o_{P}(1)$
because $\sqrt{n}\Psi_{n}(\theta_{\infty})=O_{P}(1)$ and $\mathbb{B}$ is
invertible.

The additional assumptions required for asymptotic normality (Assumption A.4)
in common choices of $L$ and $M$ (such as logistic or log-linear models) are
relatively mild. Similar to continuity, it is reasonable to assume Lipschitz
continuity and differentiability for the mapping $\theta\rightarrow
U_{\rm{stacked}}$, since $U_{\rm{stacked}}$ is a smooth function over $\Theta$
with $p_{k}(R_{k};\eta)$ bounded away from $0$ and $1$ if needed. Finiteness
of $\mathbb{B}$ is linked to integrability and finding a dominating function
for $\theta\rightarrow U_{\rm{stacked}}$, and requires that $K$ and the
variables contributing to $U_{k}$ cannot be too large. Invertibility of
$\mathbb{B}$ follows from the uniqueness of our minimizer $\theta_{\infty}$.
The main concern is the finiteness of $\mathbb{C}$ and the square-
integrability of the Lipschitz constant, which further limits the size of $K$
and the variables contributing to $U_{k}$.

###  4.4 Double Robustness

Our next task is to connect $\beta_{\infty}$, for which $\beta_{n}$ is
consistent, to the marginalized lag effect $\xi_{\rm{marg}}(k,\ell,s)$. We can
make this connection if some of our various modeling assumptions is correct.
All told, there are four different models that we wish were valid for some
$\theta\in\Theta$:

  * •

Outcome model, viz,

$\displaystyle g(R_{k})^{\prime}\alpha+A_{k}f(S_{k})^{\prime}\beta=\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}\,|\,R_{k},A_{k},T_{k+\ell}<\infty],$
(4)

  * •

Conditional probability model in the numerator of the weights, viz,

$\displaystyle q_{k}(S_{k};\xi)=\mathop{\rm
I\kern-1.99997ptE}\nolimits[A_{k}\,|\,S_{k},T_{k+\ell}<\infty],$ (5)

  * •

Conditional probability model in the denominator of the weights, viz,

$\displaystyle p_{k}(R_{k};\eta)=\mathop{\rm
I\kern-1.99997ptE}\nolimits[A_{k}\,|\,R_{k},T_{k+\ell}<\infty],$ (6)

  * •

Marginalized lag effect, viz,

$\displaystyle f(S_{k})^{\prime}\beta=\zeta_{\rm{marg}}(k,\ell,S_{k}),$ (7)

######  Theorem 3 (Double robustness)

Suppose positivity (Assumption 4) and the conditions for asymptotic normality
(Theorem 2) hold. If either the outcome model is correct (Equation 4) or the
conditional probability model in the denominator of weights (Equation 6) is
correct for $\theta_{\infty}\in\Theta$, then the asymptotic limit of our
estimator $\beta_{\infty}$ is related to the lag effect
$\zeta_{\rm{marg}}(k,\ell,S_{k})$ according to the linear equation:

$\mathbb{G}\beta_{\infty}=g$

where

$\displaystyle\mathbb{G}$ $\displaystyle=\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[q(S_{k};\xi_{\infty})(1-q(S_{k};\xi_{\infty}))f(S_{k})f(S_{k})^{\prime}\,\big{|}\,T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty)$ $\displaystyle g$
$\displaystyle=\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[q(S_{k};\xi_{\infty})(1-q(S_{k};\xi_{\infty}))\zeta_{\rm{marg}}(k,\ell,S_{k})f(S_{k})|T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty).$

If, in addition, $\mathbb{G}$ is invertible and the lag effect model is
correct (Equation 7) for some
$\theta_{*}=(\xi_{*},\eta_{*},\alpha_{*},\beta_{*})\in\Theta$, then

$\beta_{\infty}=\beta_{*}.$

The proof for this theorem can be found in Appendix C. There are several
points to consider. First, $\beta_{\infty}$ can be viewed as a weighted
average of the causal effect $\zeta_{\rm{marg}}(k,\ell,S_{k})$, even if the
lag effect model is incorrect. This implies that we have a third form of
robustness where, although we may not precisely recover our target effect, we
can still obtain a relevant causal effect. We use a weighting scheme based on
$q_{k}(S_{k};\xi_{\infty})(1-q_{k}(S_{k};\xi_{\infty}))$. This idea is not
novel, as there exists a recent method in causal inference known as overlap
weights. This method weights the average causal effect conditional on a
variable $x$, denoted by $\mu(x)$, by $e(x)(1-e(x))$, where $e(x)$ is the
propensity score conditional on $x$ (Li et al. 2019). Our weighting scheme
gives greater importance to jobs where $S_{k}$ is balanced between
intervention groups. Such jobs have intervention decisions that are
essentially random, allowing for greater comparability between the groups.
Last, double robustness is hiding an important point. The error terms in the
proof are directly related to the errors in both the conditional probability
model in the denominator (Equation 6) and the outcome model (Equation 4). Even
if both models are incorrect but are good approximations,
$\mathbb{G}\beta_{\infty}\approx g.$

###  4.5 Semiparametric Efficiency

The theorem on asymptotic normality (Theorem 2) provides some notion of
efficiency. It shows that the estimators’ variance takes the form of the
standard sandwich estimator
$\frac{1}{n}\mathbb{B}^{-1}\mathbb{C}(\mathbb{B}^{-1})^{\prime}$. This
observation enables us to estimate the variance of our estimator, along with
standard errors, by using their empirical counterpart:

$\displaystyle\mathbb{B}_{n}=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}\left[\nabla_{\theta}U_{\rm{stacked}}\big{|}_{\theta=\theta_{n}}\right];\qquad\mathbb{C}_{n}=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}\left[U_{\rm{stacked}}U_{\rm{stacked}}^{\prime}\big{|}_{\theta=\theta_{n}}\right].$

One important question is whether there is a better way to estimate the
inferential target $\beta$. This question is often asked in terms of
asymptotic efficiency, which captures the extent to which our estimator for
$\beta$ has small asymptotic variance. To compute the asymptotic variance of
our estimator of $\beta$, we can use the matrix
$\frac{1}{n}\mathbb{B}^{-1}\mathbb{C}(\mathbb{B}^{-1})^{\prime}$. In
particular, the asymptotic variance of our estimator for $\beta$ is given by
the $\beta$ portion of the formula, which can be written as:

$\frac{1}{n}\left(\nabla_{\beta}\theta\cdot\mathbb{B}^{-1}\right)\mathbb{C}\left(\nabla_{\beta}\theta\cdot\mathbb{B}^{-1}\right)^{\prime}.$

Although our estimator is not generally optimal, we offer the following
evidence to suggest that it is a sensible choice:

######  Theorem 4 (Efficient score)

Consider the special case when
$S_{k}=R_{k}=\mathop{{\rm{Pa}}}\nolimits(Y_{k+\ell})\setminus A_{k}$ and the
semi-parametric model ${\cal M}$ described in Definition A.1. Let $\mathop{\rm
I\kern-1.99997ptE}\nolimits_{0}$ and $\rm{Var}_{0}$ denote expectation and
variance with respect to the correct model in ${\cal M}.$ Then, the efficient
score function for $\beta$ with respect to ${\cal M}$ is

$\displaystyle
S_{\rm{eff}}(\beta)=\sum_{k=1}^{K-\ell}\frac{f(R_{k})}{\sigma_{k}(R_{k})}\left(Y_{k+\ell}-\mu_{k}(R_{k})-A_{k}f(R_{k})^{\prime}\beta\right)\left(A_{k}-\rho_{k}(R_{k})\right)$

where

$\displaystyle\rho_{k}(R_{k})$ $\displaystyle=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{0}\left[A_{k}|R_{k},T_{k+\ell}<\infty\right]$
$\displaystyle\mu_{k}(R_{k})$ $\displaystyle=\mathop{\rm
I\kern-1.99997ptE}\nolimits_{0}[Y_{k+\ell}|A_{k}=0,R_{k},T_{k+\ell}<\infty]$
$\displaystyle\sigma_{k}(R_{k})$
$\displaystyle=(1-\rho_{k}(R_{k})){\rm{Var}}_{0}[Y_{k+\ell}|A_{k}=1,R_{k},T_{k+\ell}<\infty]+\rho_{k}(R_{k}){\rm{Var}}_{0}[Y_{k+\ell}|A_{k}=0,R_{k},T_{k+\ell}<\infty].$

The proof of the theorem above is presented in Appendix D. It is worth noting
that the condition $S_{k}=R_{k}$ implies $W_{k}=1$, based on Equation 3. The
form of the efficient score inspired us to define the estimating equation
(Equation 1). Specifically, under the conditions stated in the theorem, our
proposed estimator satisfies the following expression:

$\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}\left[\sum_{k=1}^{K-\ell}f(R_{k})\left(Y_{k+\ell}-g(R_{k})^{\prime}\alpha-
A_{k}f(R_{k})^{\prime}\beta\right)\left(A_{k}-p_{k}(R_{k};\eta)\right)\right].$

Upon comparing this expression to the one in the theorem, we observe that
although we do not have precise knowledge of $\mu_{k}(R_{k})$, we utilize
$g(R_{k})^{\prime}\alpha$ as an approximation for $\mu_{k}(R_{k})$. Similarly,
while we may not have exact information about $\rho_{k}(R_{k})$, we use
$p_{k}(R_{k};\eta)$ as an approximation for $\rho_{k}(R_{k})$. These ideas are
then used for the general case when the conditioning set $S_{k}$ is a proper
subset of $R_{k}$ or when $R_{k}$ is a proper subset of
$\mathop{{\rm{Pa}}}\nolimits(Y_{k+\ell})\setminus A_{k}.$

##  5 Case Study

We apply our developed framework to observational data obtained from a large
academic hospital in the Midwest, where a split-flow model has been
implemented in the ED. In this increasingly adopted flow model, nurses at
triage are replaced by physicians who assess whether patients should be
directed to a vertical area for low-acuity cases or a traditional room.
Existing empirical evidence primarily focuses on evaluating the overall
effects of split-flow on ED operations and patient outcomes (Konrad et al.
2013, Wiler et al. 2016, Garrett et al. 2018). However, the specific impact of
the decision to assign a patient to the vertical area on subsequent patient
outcomes remains uncertain.

To estimate lag effects, we use data on $n=27831$ encounters with an ED split
flow model between November 1, 2016, and September 26, 2018. Since split-flow
is operated within a specific time window each day, we treated encounters on a
given day as independent panels. We excluded the 5 days that had fewer than 3
ED split-flow visits. Characteristics $X_{k}$ of patient $k$ within a panel
consisted of

  * •

Demographic information on age, gender, and race,

  * •

ED census, representing the number of patients in the ED upon arrival,

  * •

Chief complaint, categorized into four common complaints: abdominal pain,
chest pain, dyspnea, and headache. There was an additional category to capture
the remaining complaints;

  * •

Comorbidity factors including the Hierarchical Condition Category (HCC) score
(2019 Risk Adjustment Factors and Payment Rates 2018) and binary indicators of
congestive heart failure, hypertension, obesity, diabetes with and without
complications, and hypertension.

The intervention $A_{k}$ of patient $k$ is whether the patient is assigned to
the vertical area ($A_{k}=1$) or a traditional bed ($A_{k}=0$). We considered
four outcomes $Y_{k}$, with the estimation procedure repeated for each
outcome. These outcomes included time to disposition after being roomed (in
minutes), number of tests performed such as electrocardiograms and radiology
scans, admission decision, and the patient routing decision to either vertical
area or fast track area (which is simply $A_{k}$).

The lag $\ell$ was fixed at 1. We defined $R_{k}$ to include characteristics
$X_{k+\ell}$ of the future patient. To model the conditional probability
$p_{k}(R_{k};\eta)$ in the denominator of the weights, we performed GEE
regression with a logit link function and an independent working correlation
matrix. The regression model included the outcome $A_{k}$ as the dependent
variable. Independent variables included the current patient’s characteristics
$X_{k}$, the future patient’s characteristics $X_{k+\ell}$, and a fixed number
$i$ of lagged decisions $A_{k-i},\ldots,A_{k-1}$. Quadratic and cubic splines
were considered for age of the current patient with knots placed at tertiles
in the age distribution (i.e., 30, 46, and 62 years). Model comparison was
conducted using the Quasi-likelihood under the Independence model Criterion
(QICu) (Pan 2001), where different choices in $i$ and age splines were
evaluated. The model with the smallest QICu value was selected, which
consisted of $i=4$ and quadratic splines for age.

We explored various choices for $S_{k}$. We began with $S_{k}$ being an empty
set ($\emptyset$). To demonstrate the flexibility of considering different
$S_{k}$, we then set $S_{k}$ to include each of the comorbidity variables for
the future patient, which are available in $X_{k+\ell}$. To model the
conditional probability $q_{k}(S_{k};\xi)$ in the numerator of the weights, we
performed GEE regression with a logit link function and an independent working
correlation matrix. In the regression model, we treated $A_{k}$ as the
dependent variable. When $S_{k}$ was not empty, it was included as an
independent variable in the model.

For each outcome and each $S_{k}$, we estimated the inferential target of
$\beta$. In the context of Theorem 3, $f(S_{k})^{\prime}\beta$ corresponds to
the marginalized lag effect $\zeta_{\rm{marg}}(k,\ell,S_{k})$. To achieve
this, we solved the estimating equation presented in Equation (1) for a given
baseline model $g(R_{k})^{\prime}\alpha$ and marginalized lag effect model
$f(S_{k})^{\prime}\beta$. The term $f(S_{k})$ included a constant term of $1$
to represent the main intervention effect, and when $S_{k}$ was one of the
comorbidity variables, it incorporated $S_{k}$ as well. Our baseline model
included a constant term of $1$ to represent the intercept, as well as the
characteristics of the current ($X_{k}$) and future ($X_{k+\ell}$) patient.
Additionally, the requisite term $f(S_{k})q_{k}(S_{k};\xi_{n})$ was included,
where $\xi_{n}$ represents the estimate obtained from performing GEE
regression to model $q_{k}(S_{k};\xi)$. As mentioned earlier, we tackled the
estimating equation by employing an equivalent procedure of GEE regression
with an identity link and independent working correlation matrix.

Tables 1 and 2 reports the main effect $\beta$ when $S_{k}=\emptyset$ and the
interaction term $f(1)^{\prime}\beta-f(0)^{\prime}\beta$ when
$S_{k}\neq\emptyset$. Standard errors were calculated using the geex package
(Saul and Hudgens 2020). The implementation details are in Appendix E.
Estimates and standard errors were used to construct Wald 95% confidence
intervals and perform Wald hypothesis tests of the null hypothesis that the
true value of the estimate is zero. Significance was considered $P<.05$. Given
the exploratory nature of this investigation, marginal significance was
considered $P<.10$.

Table 1: Estimates for main effect ($S_{k}=\emptyset$) or interaction term ($S_{k}\neq\emptyset$) of sending the current patient to the vertical area on the time to disposition and number of tests for the next patient. We report $P$ values for Wald tests. | Time to disposition (min.) | Number of tests
---|---|---
Variable $S_{k}$ | Estimate (95% CI) | $P$ | Estimate (95% CI) | $P$
$\emptyset$ | 3.6 (-1.1, 8.3) | .11 | 0.04 (0.01, 0.09) | .02
Congestive heart failure | -3.1 (-22.3, 16.1) | .74 | -0.06 (-0.2, 0.1) | .56
Hypertension | -1.9 (-11.0, 7.2) | .68 | -0.03 (-0.1, 0.05) | .43
Obesity | 2.9 (-9.4, 15.4) | .63 | 0.02 (-0.1, 0.2) | .70
Diabetes with complications | -10.5 (-25.4, 4.4) | .16 | -0.01 (-0.1, 0.1) | .88
Diabetes without complications | -7.3 (-22.2, 7.5) | .33 | 0.06 (-0.1, 0.2) | .45
HCC | -1.9 (-4.1, 0.3) | .09 | 0.01 (-0.01, 0.03) | .47
Table 2: Estimates for main effect ($S_{k}=\emptyset$) or interaction term ($S_{k}\neq\emptyset$) of sending the current patient to the vertical area on the admission decision and vertical area decision for the next patient. We report $P$ values for Wald tests. | Admission decision (%) | Vertical area decision (%)
---|---|---
Variable $S_{k}$ | Interaction (95% CI) | $P$ | Interaction (95% CI) | $P$
$\emptyset$ | 0.9 (-0.8, 2.7) | .29 | 2.0 (0.3, 3.6) | .01
Congestive heart failure | 1.6 (-5.5, 8.8) | .65 | -0.7 (-5.5, 4.1) | .76
Hypertension | -3.7 (-7.4, 0.05) | .05 | 0.3 (-2.6, 3.3) | .81
Obesity | -0.5 (-5.0, 3.9) | .81 | 0.1 (-4.1, 4.4) | .94
Diabetes with complications | 1.0 (-5.2, 7.3) | .75 | -0.1 (-5.2, 4.8) | .95
Diabetes without complications | -4.0 (-11.0, 2.9) | .25 | 0.6 (-5.2, 6.5) | .83
HCC | -0.1 (-1.1, 0.7) | .74 | -0.3 (-1.0, 0.4) | .46

Our results suggest sequential bias in the decision to send a patient to a
vertical area over a traditional bed. In Table 1, we find that sending the
current patient to vertical area leads to a significant increase of $0.04$
(95% CI: $[0.01,0.09]$) tests on average for the next patient. Regarding time
to disposition, our findings indicate a marginally significant interaction
between the next patient’s HCC score and the decision to send the current
patient to the vertical area. Although, in general, sending the current
patient to vertical area leads to an estimated non-significant increase in
average time to disposition of 3.6 minutes for the subsequent patient, this
increase is 1.9 (95% CI: $[-0.3,4.1]$) minutes shorter when the subsequent
patient has an HCC score of 1 compared to an HCC score of 0. In Table 2, we
find that sending the current patient to vertical area leads to a significant
increase of 2% (95% CI: $[0.3\%,3.6\%]$) in the next patient’s probability of
being assigned to vertical area. Further, we observe a marginally significant
interaction between the next patient’s hypertension status and the decision to
send the current patient to the vertical area. Although, in general, sending
the current patient to vertical area lead to an estimated non-significant
increase of 0.9% in the next patient’s admission rate, this rate increase is
estimated to be -3.7% (95% CI: $[-7.4\%,0.05\%]$) lower when the patient has
hypertension compared to when they do not have hypertension.

##  6 Conclusion

We introduced a new causal inference framework to examine sequential bias in
stochastic service systems from observational data. Sequential bias has been
widely observed across diverse domains, including sporting competitions, stock
markets, and judging decisions. If we can accurately measure this bias, we can
pave the way for enhancing decision-making in stochastic service systems
through stochastic modeling and optimization. With this goal in mind, we had
to address two obstacles to performing causal inference: interference and a
random number of jobs. To handle interference, we drew upon insights from
DTRs, building specifically upon the work in Robins (1997), Boruvka et al.
(2018), Qian et al. (2021). Handling the random number of jobs introduced
challenges across modeling, identification, and estimation. These challenges
include building a causal model that permits a random number of variables,
identifying the resulting causal effects, and deriving estimation properties
such as double robustness and an efficient score. Our innovative solutions to
these challenges represent the core contributions of our work, pushing the
boundaries of causal inference forward.

To tackle the issue of modeling a random number of jobs, we conceptualized the
decision-making process as a MPP and integrated it within the well-established
causal inference framework established by Richardson and Robins (2013) and
Malinsky et al. (2019). By employing the MPP, we were able to effectively
account for the number of jobs. However, this required working with an
infinite number of random variables, posing certain technical difficulties.
Ensuring the integrability of our proposed estimating equations and justifying
the exchange between differentiation and an infinite sum in our proof of
semiparametric efficiency were among the difficulties we encountered and
successfully navigated. These successes not only help us quantify sequential
bias, but also offer a way to effectively handle a random number of decisions
in a DTR.

Based on the MPP, we introduced two causal effects as measures of sequential
bias. These effects differ from their counterpart in the DTR literature
(Boruvka et al. 2018, Qian et al. 2021) for their inclusion of variables
collected after the intervention in the conditioning set. This inclusion is
unavoidable because, in the presence of a random number of jobs, measuring the
causal impact that intervening on a current job has on a future job requires
that the future job exists. In other words, we always condition on a specific
post-intervention variable, namely the binary indicator of whether the future
job exists. Moreover, the inclusion of post-intervention variables brings an
additional advantage of modeling the future job’s outcome using its own
characteristics, rather than the characteristics of the current or past job.
This approach holds promise for achieving more efficient estimators. Indeed,
we derived an efficient score for those specific lag effects that conditioned
on _every_ post-intervention variable that is a parent of the future job’s
outcome. Without conditioning on post-intervention variables, our proof would
only work for effects on the outcome of the currently intervened job, instead
of a future job. This consideration may explain why Qian et al. (2021), which
did not condition on post-intervention variables, provided an efficient score
only in the special case of their causal effect when the outcome was the most
immediate.

Nevertheless, it is crucial to acknowledge a caveat when conditioning on post-
intervention variables: the interpretation of the lag effect becomes more
nuanced. To address the confusion surrounding conditioning on post-
intervention variables, Pearl (2015) extensively discusses this topic in a
dedicated paper. He highlights the importance of considering whether to
condition on the observed value of post-intervention variables or the
hypothetical value that these variables would achieve under the contrasting
scenarios. These considerations lead to different research questions. In our
study, the question was clear: to identify the lag effect, we needed to
condition on the hypothetical value rather than the observed value.
Consequently, the resulting lag effect is interpreted as an average contrast
in a future job’s outcome between intervening and not intervening on a current
job among those events whose post-treatment variables would achieve a certain
value regardless of the intervention. For instance, we may compare the future
job’s outcome when intervening or not intervening on a current job, focusing
on events where the future job exists in either case.

Applying our framework to specific settings comes with several limitations to
consider. One such limitation is that, like many other causal inference
approaches, we rely on assumptions that cannot be directly verified from data.
These assumptions involve considering multiple hypothetical scenarios, of
which only one is actually realized. Therefore, the validity of our approach
in a given application depends on how well the application can be accurately
represented using our proposed causal model. A specific issue to be mindful
of, because it could compromise the model’s validity, is an unobserved common
cause of the future job’s outcome and the current job’s outcome. Thus, careful
consideration of such confounding factors is necessary to ensure the validity
of our approach.

Another limitation is the unknown assignment mechanism for the intervention,
which needs to be estimated from the available data. Since our estimator is
doubly-robust, this implies that if both the probability model of assignment
and the outcome model are incorrect, there is a risk of biased estimates. A
related issue is the occurrence of intervention assignments with extremely low
probabilities, which we then divide by during estimation. The accuracy of
these probabilities can impact on estimation. Furthermore, such probabilities
suggest that certain jobs are very unlikely to experience one of the
intervention conditions, which violates our assumption of positivity. One
proposed solution is to exclude cases with extremely low probabilities (Lee et
al. 2011). Alternatively, overlap weights can be used to smoothly reduce the
influence of observations with extreme probabilities (Li et al. 2019).
Encouragingly, our investigation of double robustness demonstrates that by
adding a term in the numerator of our weights related to the probability of
assignment, our estimator effectively produces a weighted causal effect
similar to one that used overlap weights.

A final limitation is the assumption of non-interference among panels, i.e.
that the outcomes within a panel are solely determined by the decisions in
that panel, independent of decisions made in other panels. However, this
assumption may not hold if the decision-maker for one panel of data were to
coordinate their decisions with the decision-maker on another panel. For
example, physicians in the ED can potentially impact each other’s decisions,
making it inappropriate to consider patients treated by different doctors
during the same time period as arising from distinct panels.

To conclude our investigation, we estimated lag effects for a specific
decision-making scenario in the ED, offering a concrete illustration of how
our proposed causal inference framework can be applied. We provide step-by-
step guidance on how to implement our method using commonly available
statistical software. The estimation of standard errors may be the most
challenging aspect, and we include code snippets in the Appendix to facilitate
this process. Our analysis reveals a noteworthy finding regarding the
influence of a current patient’s routing decision on the next patient.
Specifically, our data analysis demonstrates that routing the current patient
to a vertical area significantly increases both the number of tests performed
and the probability of the next patient being routed to the vertical area.
These findings provide compelling support for the existence of sequential bias
in a setting of paramount importance and profound consequences.

## References

  * (1)
  * 2019 Risk Adjustment Factors and Payment Rates (2018) 2019 Risk Adjustment Factors and Payment Rates (2018), https://www.cms.gov/Medicare/Health-Plans/MedicareAdvtgSpecRateStats/Downloads/RTC-Dec2018.pdf.  Accessed: March 23, 2023.
  * Alagoz et al. (2010) Alagoz, O., Hsu, H., Schaefer, A. J. and Roberts, M. S. (2010), ‘Markov decision processes: a tool for sequential decision making under uncertainty’, Medical Decision Making 30(4), 474–483.
  * Aronow and Samii (2017) Aronow, P. M. and Samii, C. (2017), ‘Estimating average causal effects under general interference, with application to a social network experiment’, The Annals of Applied Statistics 11(4).
  * Boruvka et al. (2018) Boruvka, A., Almirall, D., Witkiewitz, K. and Murphy, S. A. (2018), ‘Assessing Time-Varying Causal Effect Moderation in Mobile Health’, Journal of the American Statistical Association 113(523), 1112–1121.
  * Chakraborty and Moodie (2013) Chakraborty, B. and Moodie, E. E. (2013), ‘Statistical methods for dynamic treatment regimes’, Springer-Verlag. doi 10, 978–1.
  * Chakraborty and Murphy (2014) Chakraborty, B. and Murphy, S. A. (2014), ‘Dynamic Treatment Regimes’, Annual Review of Statistics and Its Application 1(1), 447–464.
  * Chen et al. (2016) Chen, D. L., Moskowitz, T. J. and Shue, K. (2016), ‘Decision Making Under the Gambler’s Fallacy: Evidence from Asylum Judges, Loan Officers, and Baseball Umpires*’, The Quarterly Journal of Economics 131(3), 1181–1242.
  * Cochran et al. (2023) Cochran, A., Maronge, J. M., Victory, A., Hoel, S., McInnis, M. G., Thomas, E. B. et al. (2023), ‘Mobile acceptance and commitment therapy in bipolar disorder: Microrandomized trial’, JMIR Mental Health 10(1), e43164.
  * Cox (1958) Cox, D. R. (1958), Planning of experiments, Planning of experiments, Wiley, Oxford, England.  Pages: 308.
  * David et al. (2022) David, J. C., Cochran, A. L., Patterson, B. W. and Zayas-Caban, G. (2022), ‘Evaluation of a split flow model for the emergency department’, arXiv preprint arXiv:2202.00736 .
  * Didelez (2008) Didelez, V. (2008), ‘Graphical models for marked point processes based on local independence’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70(1), 245–264.
  * Didelez (2015) Didelez, V. (2015), Causal reasoning for events in continuous time: A decision-theoretic approach., in ‘ACI@ UAI’, pp. 40–45.
  * Gao et al. (2021) Gao, T., Subramanian, D., Bhattacharjya, D., Shou, X., Mattei, N. and Bennett, K. P. (2021), ‘Causal inference for event pairs in multivariate point processes’, Advances in Neural Information Processing Systems 34, 17311–17324.
  * Garrett et al. (2018) Garrett, J. S., Berry, C., Wong, H., Qin, H. and Kline, J. A. (2018), ‘The effect of vertical split-flow patient management on emergency department throughput and efficiency’, The American Journal of Emergency Medicine .
  * Gill and Robins (2001) Gill, R. D. and Robins, J. M. (2001), ‘Causal Inference for Complex Longitudinal Data: The Continuous Case’, The Annals of Statistics 29(6), 1785–1811.  Publisher: Institute of Mathematical Statistics.
  * Goldbach et al. (2022) Goldbach, C., Sickmann, J. and Pitz, T. (2022), ‘Sequential decision bias – evidence from grading exams’, Applied Economics 54(32), 3727–3739.
  * Guo et al. (2021) Guo, F. R., Richardson, T. S. and Robins, J. M. (2021), ‘Discussion of ‘estimating time-varying causal excursion effects in mobile health with binary outcomes”, Biometrika 108(3), 541–550.
  * Hinderer et al. (2016) Hinderer, K., Rieder, U. and Stieglitz, M. (2016), Dynamic optimization, Springer.
  * Hong and Raudenbush (2006) Hong, G. and Raudenbush, S. W. (2006), ‘Evaluating Kindergarten Retention Policy: A Case Study of Causal Inference for Multilevel Observational Data’, Journal of the American Statistical Association 101(475), 901–910.
  * Horvitz and Thompson (1952) Horvitz, D. G. and Thompson, D. J. (1952), ‘A generalization of sampling without replacement from a finite universe’, Journal of The American Statistical Association 47(260), 663–685.
  * Hudgens and Halloran (2008) Hudgens, M. G. and Halloran, M. E. (2008), ‘Toward Causal Inference With Interference’, Journal of the American Statistical Association 103(482), 832–842.
  * Jacobsen and Gani (2006) Jacobsen, M. and Gani, J. (2006), ‘Point process theory and applications: marked point and piecewise deterministic processes’.
  * Klasnja et al. (2015) Klasnja, P., Hekler, E. B., Shiffman, S., Boruvka, A., Almirall, D., Tewari, A. and Murphy, S. A. (2015), ‘Micro-Randomized Trials: An Experimental Design for Developing Just-in-Time Adaptive Interventions’, Health Psychology : Official Journal of the Division of Health Psychology, American Psychological Association 34(0), 1220–1228.
  * Konrad et al. (2013) Konrad, R., DeSotto, K., Grocela, A., McAuley, P., Wang, J., Lyons, J. and Bruin, M. (2013), ‘Modeling the impact of changing patient flow processes in an emergency department: Insights from a computer simulation study’, Operations Research for Health Care 2(4), 66–74.
  * 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(1), 1225.
  * Lavori and Dawson (2000) Lavori, P. W. and Dawson, R. (2000), ‘A design for testing clinical strategies: Biased adaptive within-subject randomization’, Journal of the Royal Statistical Society. Series A (Statistics in Society) 163(1), 29–38.
  * Lavori and Dawson (2004) Lavori, P. W. and Dawson, R. (2004), ‘Dynamic treatment regimes: practical design considerations’, Clinical trials 1(1), 9–20.
  * Lee et al. (2011) Lee, B. K., Lessler, J. and Stuart, E. A. (2011), ‘Weight trimming and propensity score weighting’, PloS one 6(3), e18174.
  * Li et al. (2019) Li, F., Thomas, L. E. and Li, F. (2019), ‘Addressing extreme propensity scores via the overlap weights’, American Journal of Epidemiology 188(1), 250–257.
  * Li et al. (2023) Li, Z., Chen, J., Laber, E., Liu, F. and Baumgartner, R. (2023), ‘Optimal treatment regimes: a review and empirical comparison’, International Statistical Review .
  * Malinsky et al. (2019) Malinsky, D., Shpitser, I. and Richardson, T. (2019), A potential outcomes calculus for identifying conditional path-specific effects, in ‘The 22nd International Conference on Artificial Intelligence and Statistics’, PMLR, pp. 3080–3088.
  * Murphy (2003) Murphy, S. A. (2003), ‘Optimal dynamic treatment regimes’, Journal of the Royal Statistical Society. Series B (Statistical Methodology) 65(2), 331–366.
  * Murphy (2005) Murphy, S. A. (2005), ‘An experimental design for the development of adaptive treatment strategies’, Statistics in Medicine 24(19), 3179–3194.
  * Nahum-Shani et al. (2017) Nahum-Shani, I., Smith, S. N., Spring, B. J., Collins, L. M., Witkiewitz, K., Tewari, A. and Murphy, S. A. (2017), ‘Just-in-Time Adaptive Interventions (JITAIs) in Mobile Health: Key Components and Design Principles for Ongoing Health Behavior Support’, Annals of Behavioral Medicine: A Publication of the Society of Behavioral Medicine 52(6), 446–462.
  * Neyman (1923) Neyman, J. (1923), ‘On the application of probability theory to agricultural experiments. essay on principles’, Ann. Agricultural Sciences pp. 1–51.
  * Pan (2001) Pan, W. (2001), ‘Akaike’s information criterion in generalized estimating equations’, Biometrics 57(1), 120–125.
  * Pearl (2000) Pearl, J. (2000), Causality: models, reasoning, and inference, Cambridge University Press, Cambridge, U.K. ; New York.
  * Pearl (2015) Pearl, J. (2015), ‘Conditioning on post-treatment variables’, Journal of Causal Inference 3(1), 131–137.
  * Qian et al. (2021) Qian, T., Yoo, H., Klasnja, P., Almirall, D. and Murphy, S. A. (2021), ‘Estimating time-varying causal excursion effects in mobile health with binary outcomes’, Biometrika 108(3), 507–527.
  * Richardson and Robins (2013) Richardson, T. S. and Robins, J. M. (2013), ‘Single world intervention graphs (swigs): A unification of the counterfactual and graphical approaches to causality’, Center for the Statistics and the Social Sciences, University of Washington Series. Working Paper 128(30), 2013.
  * Robins (1986) Robins, J. (1986), ‘A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect’, Mathematical modelling 7(9-12), 1393–1512.
  * Robins (1994) Robins, J. M. (1994), ‘Correcting for non-compliance in randomized trials using structural nested mean models’, Communications in Statistics - Theory and Methods 23(8), 2379–2412.
  * Robins (1997) Robins, J. M. (1997), Causal Inference from Complex Longitudinal Data, in M. Berkane, ed., ‘Latent Variable Modeling and Applications to Causality’, Lecture Notes in Statistics, Springer, New York, NY, pp. 69–117.
  * Røysland (2011) Røysland, K. (2011), ‘A martingale approach to continuous-time marginal structural models’, Bernoulli 17(3), 895–915.
  * Røysland (2012) Røysland, K. (2012), ‘Counterfactual analyses with graphical models based on local independence’, The Annals of Statistics pp. 2162–2194.
  * Rubin (1974) Rubin, D. B. (1974), ‘Estimating causal effects of treatments in randomized and nonrandomized studies.’, Journal of Educational Psychology 66(5), 688.
  * Rubin (1980) Rubin, D. B. (1980), ‘Randomization analysis of experimental data: The fisher randomization test comment’, Journal of the American Statistical Association 75(371), 591–593.
  * Rubin (2005) Rubin, D. B. (2005), ‘Causal Inference Using Potential Outcomes’, Journal of the American Statistical Association 100(469), 322–331.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1198/016214504000001880.
  * Saul and Hudgens (2020) Saul, B. C. and Hudgens, M. G. (2020), ‘The calculus of m-estimation in R with geex’, Journal of Statistical Software 92(2), 1–15.
  * Shi et al. (2022) Shi, J., Wu, Z. and Dempsey, W. (2022), ‘Estimating time-varying direct and indirect causal excursion effects with longitudinal binary outcomes’, arXiv preprint arXiv:2212.01472 .
  * Sobel (2006) Sobel, M. E. (2006), ‘What Do Randomized Studies of Housing Mobility Demonstrate?: Causal Inference in the Face of Interference’, Journal of the American Statistical Association 101(476), 1398–1407.
  * Sofrygin and van der Laan (2017) Sofrygin, O. and van der Laan, M. J. (2017), ‘Semi-Parametric Estimation and Inference for the Mean Outcome of the Single Time-Point Intervention in a Causally Connected Population’, Journal of Causal Inference 5(1), 20160003.
  * Steimle and Denton (2017) Steimle, L. N. and Denton, B. T. (2017), ‘Markov decision processes for screening and treatment of chronic diseases’, Markov Decision Processes in Practice pp. 189–222.
  * Tauchen (1985) Tauchen, G. (1985), ‘Diagnostic testing and evaluation of maximum likelihood models’, Journal of Econometrics 30(1-2), 415–443.
  * Tchetgen Tchetgen et al. (2021) Tchetgen Tchetgen, E. J., Fulcher, I. R. and Shpitser, I. (2021), ‘Auto-G-Computation of Causal Effects on a Network’, Journal of the American Statistical Association 116(534), 833–844.
  * Thomas et al. (2023) Thomas, E. B. K., Sagorac Gruichich, T., Maronge, J. M., Hoel, S., Victory, A., Stowe, Z. N. and Cochran, A. (2023), ‘Mobile acceptance and commitment therapy with distressed first-generation college students: Microrandomized trial’, JMIR Mental Health 10, e43065.
  * Tsiatis (2006) Tsiatis, A. A. (2006), Semiparametric theory and missing data, Springer series in statistics, Springer, New York.
  * Tsiatis et al. (2019) Tsiatis, A. A., Davidian, M., Holloway, S. T. and Laber, E. B. (2019), Dynamic treatment regimes: Statistical methods for precision medicine, CRC press.
  * Van der Vaart (2000) Van der Vaart, A. W. (2000), Asymptotic statistics, Vol. 3, Cambridge university press.
  * Vansteelandt and Sjolander (2016) Vansteelandt, S. and Sjolander, A. (2016), ‘Revisiting g-estimation of the effect of a time-varying exposure subject to time-varying confounding’, Epidemiologic Methods 5(1), 37–56.
  * Verbitsky-Savitz and Raudenbush (2012) Verbitsky-Savitz, N. and Raudenbush, S. W. (2012), ‘Causal Inference Under Interference in Spatial Settings: A Case Study Evaluating Community Policing Program in Chicago’, Epidemiologic Methods 1(1).
  * Wang and Yin (2020) Wang, X. and Yin, L. (2020), ‘New $G$-formula for the sequential causal effect and blip effect of treatment in sequential causal inference’, The Annals of Statistics 48(1).
  * Wiler et al. (2016) Wiler, J. L., Ozkaynak, M., Bookman, K., Koehler, A., Leeret, R., Chua-Tuan, J., Ginde, A. A. and Zane, R. (2016), ‘Implementation of a front-end split-flow model to promote performance in an urban academic emergency department’, The Joint Commission Journal on Quality and Patient Safety 42(6), 271–AP4.
  * Yu and Cohen (2008) Yu, A. J. and Cohen, J. D. (2008), ‘Sequential effects: Superstition or rational behavior?’, Advances in Neural Information Processing Systems 21, 1873–1880.
  * Zayas-Cabán et al. (2016) Zayas-Cabán, G., Xie, J., Green, L. V. and Lewis, M. E. (2016), ‘Dynamic control of a tandem system with abandonments’, Queueing Systems 84, 279–293.
  * Zhang et al. (2022) Zhang, C., Mohan, K. and Pearl, J. (2022), ‘Causal Inference with Non-IID Data using Linear Graphical Models’, p. 12.

Appendix

##  Appendix A Regularity Conditions

To construct the distribution for our MPP, we make use of the Ionescu Tulcea
theorem. Letting $\oplus_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}{\cal
E}_{k}$ be the product $\sigma$-algebra in $\prod_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}\Omega_{k}$ given a sequence $(\Omega_{k},{\cal
E}_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$ of measurable spaces,
the statement of the theorem is the following:

######  Theorem A.1 (Ionescu Tulcea)

For a sequence of measurable spaces $(\Omega_{k},{\cal
E}_{k})_{k\in\mathop{\rm I\kern-1.39998ptN}\nolimits}$, a probability measure
$\mathop{\rm I\kern-1.99997ptP}\nolimits_{0}$ on $\Omega_{1}$, and Markov
kernels $\mathop{\rm I\kern-1.99997ptP}\nolimits_{k}$ from
$\prod_{i=1}^{k}\Omega_{k}$ to $\Omega_{k+1}$ for $k\in\mathop{\rm
I\kern-1.99997ptN}\nolimits$, there exists random variables $Z_{k}$ taking
values in $\Omega_{k}$ with unique probability distribution $\mathop{\rm
I\kern-1.99997ptP}\nolimits$ on $(\prod_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}\Omega_{k},\otimes_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}{\cal E}_{k})$ satisfying the following equation
for all $k\in\mathop{\rm I\kern-1.99997ptN}\nolimits$ and $B_{i}\in{\cal
E}_{i}$:

$\mathop{\rm I\kern-1.99997ptP}\nolimits(\\{Z_{i}\in
B_{i}\\}_{i=1}^{k})=\int_{B_{1}}\mathop{\rm
I\kern-1.99997ptP}\nolimits_{0}(d\omega_{1})\int_{B_{2}}\mathop{\rm
I\kern-1.99997ptP}\nolimits_{1}(d\omega_{2}|\omega_{1})\ldots\int_{B_{k}}\mathop{\rm
I\kern-1.99997ptP}\nolimits_{k}(d\omega_{k}|\omega_{1},\ldots,\omega_{k-1}).$

In order for the random variables $(T_{k},X_{k},A_{k},Y_{k})_{k\in\mathop{\rm
I\kern-1.39998ptN}\nolimits}$ to form an MPP, additional conditions on the
Markov kernels must be satisfied. These conditions ensure that the times are
strictly increasing when they are finite, and increasing otherwise, with the
marks taking on the irrelevant mark when the corresponding time is infinite.
Specifically, the following assumptions are made:

###### Assumption A.1

(Marked point process conditions) The Markov kernels from Assumption 2
satisfy:

$\displaystyle\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(T_{k})}\left((t_{k-1},\infty]\,\big{|}\,t_{1:k-1},h_{k-1}\right)$
$\displaystyle=1$ $\displaystyle\text{if }\,\,t_{k-1}<\infty$
$\displaystyle\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(T_{k})}\left(\\{\infty\\}\,\big{|}\,t_{1:k-1},h_{k-1}\right)$
$\displaystyle=1$ $\displaystyle\text{if }\,\,t_{k-1}=\infty$
$\displaystyle\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(X_{k})}\left(\mathbb{R}^{l}\,\big{|}\,t_{1:k},h_{k-1}\right)$
$\displaystyle=1$ $\displaystyle\text{if }\,\,t_{k}<\infty$
$\displaystyle\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(X_{k})}\left(\\{\Delta\\}\,\big{|}\,t_{1:k},h_{k-1}\right)$
$\displaystyle=1$ $\displaystyle\text{if }\,\,t_{k}=\infty.$

and

$\displaystyle\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(A_{k})}\left(\\{0,1\\}\,\big{|}\,x_{k},h_{k-1}\right)$
$\displaystyle=\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(Y_{k})}\left(\mathbb{R}\,\big{|}\,x_{k},a_{k},h_{k-1}\right)=1$
$\displaystyle\text{if }\,\,x_{k}\neq\Delta$ $\displaystyle\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(A_{k})}\left(\\{\Delta\\}\,\big{|}\,x_{k},h_{k-1}\right)$
$\displaystyle=\mathop{\rm
I\kern-1.99997ptP}\nolimits^{(Y_{k})}\left(\\{\Delta\\}\,\big{|}\,x_{k},a_{k},h_{k-1}\right)=1$
$\displaystyle\text{if }\,\,x_{k}=\Delta.$

For estimation, we need several assumptions for the proposed estimator to be
integrable (Lemma 2), consistent (Theorem 1), and asymptotically normal
(Theorem 2). The next assumption ensures that $U_{\rm{stacked}}$ is
integrable.

###### Assumption A.2

(Integrability conditions) The random number of jobs has finite expectation,
viz, $\mathop{\rm I\kern-1.99997ptE}\nolimits[K]<\infty$. We also have the
following bound:

$\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\begin{Vmatrix}\begin{bmatrix}U_{k}\\\
L_{k}\\\ M_{k}\end{bmatrix}\end{Vmatrix}\mathbbm{1}_{\\{k\leq K\\}}\right]\leq
C\mathop{\rm I\kern-1.99997ptP}\nolimits(k\leq K)$

for some constant $C$.

In addition to assumptions above, we also invoke the assumptions below for
consistency and asymptotic normality. The following assumptions for
consistency incorporate the assumptions in Lemma 1 for Tauchen (1985), with
the exception that the separability condition is replaced with a more direct
condition on measurability of a function involving a supremum:

###### Assumption A.3

(Consistency conditions) In addition to the conditions needed for
integrability (Assumption A.2), the set $\Theta$ is compact, $\mathop{\rm
I\kern-1.99997ptE}\nolimits[U_{\rm{stacked}}]$ has a unique root
$\theta_{\infty}$ in $\Theta$, and almost surely $\mathop{\rm
I\kern-1.99997ptE}\nolimits_{n}[U_{\rm{stacked}}]$ has a root $\theta_{n}$ in
$\Theta$. Take $Z:=(X_{k},A_{k},Y_{k})_{k\in\\{1,\ldots,K\\}}$. Let $\psi$ be
the mapping from $\theta$ and $Z$ to $\psi(Z,\theta)=U_{\rm{stacked}}$ and
assume $\psi$ is

  * •

almost surely continuous for each $\theta\in\Theta$

  * •

almost surely bounded above by a measurable function $d(z)$ that is invariant
to $\theta$ and has finite mean

  * •

measurable for each $\theta\in\Theta$

  * •

is defined such that

$\sup_{\theta\in\Theta:\begin{Vmatrix}\theta-\theta^{*}\end{Vmatrix}<\delta}\begin{Vmatrix}U_{\rm{stacked}}-U_{\rm{stacked}}\big{|}_{\theta=\theta^{*}}\end{Vmatrix}$

is measurable for all $\delta>0$ and $\theta^{*}\in\Theta$.

The hardest step in the proof of consistency (Theorem 1) is showing a uniform
law of large numbers from these conditions, for which we refer the reader to
the proof of Lemma 1 in Tauchen (1985). Asymptotic normality (Theorem 2)
requires yet stronger conditions, which are given as followed:

###### Assumption A.4

(Asymptotic normality conditions) In addition to the conditions needed for
integrability (Assumption A.2) and consistency (Assumption A.3), we have

  * •

Measurable $\dot{\psi}(z)$ exists with $\mathop{\rm
I\kern-1.99997ptE}\nolimits[\dot{\psi}(Z)^{2}]<\infty$ such that for
$\theta_{1},\theta_{2}$ in $\Theta$ and in a neighborhood of $\theta_{\infty}$

$\begin{Vmatrix}\psi(z,\theta_{1})-\psi(z,\theta_{2})\end{Vmatrix}\leq\dot{\psi}(x)\begin{Vmatrix}\theta_{1}-\theta_{2}\end{Vmatrix}$

  * •

The matrices

$\displaystyle\mathbb{B}=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\nabla_{\theta}U_{\rm{stacked}}\big{|}_{\theta=\theta_{\infty}}\right];\qquad\mathbb{C}=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[U_{\rm{stacked}}U_{\rm{stacked}}^{\prime}\big{|}_{\theta=\theta_{\infty}}\right],$

exist and are invertible.

##  Appendix B Proof of Integrability (Lemma 2)

With the assumptions for integrability stated in the last section, we are
ready to provide the proof of integrability (Lemma 2):

#### Proof.

Our proof is for $U$. The same argument holds for $L$ and $M$. Note that the
triangle inequality implies that
$\begin{Vmatrix}U\end{Vmatrix}\leq\sum_{k=1}^{K}\begin{Vmatrix}U_{k}\end{Vmatrix}$
almost surely. As a result, we have that

$\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\begin{Vmatrix}U\end{Vmatrix}\right]=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\begin{Vmatrix}\sum_{k=1}^{K}U_{k}\end{Vmatrix}\right]\leq\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\sum_{k=1}^{K}\begin{Vmatrix}U_{k}\end{Vmatrix}\right].$

Next note that we can write $\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\sum_{k=1}^{K}\begin{Vmatrix}U_{k}\end{Vmatrix}\right]$
as $\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\sum_{k=1}^{\infty}\begin{Vmatrix}U_{k}\end{Vmatrix}\mathbbm{1}_{\\{k\leq
K\\}}\right]$. Applying the monotone convergence theorem to this last
expression allows us to write $\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\sum_{k=1}^{\infty}\begin{Vmatrix}U_{k}\end{Vmatrix}\mathbbm{1}_{\\{k\leq
K\\}}\right]=\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\begin{Vmatrix}U_{k}\end{Vmatrix}\mathbbm{1}_{\\{k\leq
K\\}}\right]$. The second part of Assumption A.2 is that

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\begin{Vmatrix}U_{k}\end{Vmatrix}\mathbbm{1}_{\\{k\leq
K\\}}\right]\leq C\mathop{\rm I\kern-1.99997ptP}\nolimits(k\leq K),$

which implies

$\displaystyle\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\begin{Vmatrix}U_{k}\end{Vmatrix}\mathbbm{1}_{\\{k\leq
K\\}}\right]\leq C\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptP}\nolimits(k\leq K)=C\mathop{\rm
I\kern-1.99997ptE}\nolimits[K].$

The first part of Assumption A.2 gives the finiteness of $\mathop{\rm
I\kern-1.99997ptE}\nolimits[K]$, which completes the proof.

##  Appendix C Proof of Double Robustness (Theorem 3)

Here, we provide the proof for Theorem 3, which relates our estimator to
$\beta_{\infty}$ to our target effect $\zeta_{\rm{marg}}(k,\ell,S_{k})$ in
situations when some of our models are correct.

#### Proof.

The proof of this theorem involves understanding what it means for
$\theta_{\infty}$ to be a solution for $\mathop{\rm
I\kern-1.99997ptE}\nolimits[U_{\rm{stacked}}].$ We accomplish this in several
steps.

Step 1: Re-write the estimating equation. We first note that

$\mathop{\rm I\kern-1.99997ptE}\nolimits[U_{\rm{stacked}}]=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\sum_{k=1}^{K-\ell}U_{k}\right]=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[\sum_{k=1}^{\infty}U_{k}\mathbbm{1}_{\\{k\leq
K\\}}\right]=\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[U_{k}|T_{k}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty).$

We were able to swap the infinite sum with the expectation by utilizing the
conditions for integrability (Assumption A.2) and the dominated convergence
theorem. In particular, we established the finiteness of the expectation of
$\begin{Vmatrix}U_{\rm{stacked}}\end{Vmatrix}$ in Lemma 2, which dominates
$U_{\rm{stacked}}.$ Further, the solution $\theta_{\infty}$ to the estimating
equation:

$\mathop{\rm
I\kern-1.99997ptE}\nolimits[U_{\rm{stacked}}]=\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[W_{k}\left(Y_{k+\ell}-g(R_{k})^{\prime}\alpha-
A_{k}f(R_{k})^{\prime}\beta\right)\begin{bmatrix}g(R_{k})\\\
A_{k}f(S_{k})\end{bmatrix}\bigg{|}T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty)=0$

is also a solution to

$\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[W_{k}\left(Y_{k+\ell}-g(R_{k})^{\prime}\alpha-
A_{k}f(S_{k})^{\prime}\beta\right)(A_{k}-q_{k}(S_{k};\xi))f(S_{k})\,\big{|}\,T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty)=0.$

This can be shown by subtracting the rows of $U_{k}$ involving
$q_{k}(S_{k};\xi)f(S_{k})$ in $g(R_{k})$ from the rows involving
$A_{k}f(S_{k}),$ which is a direct consequence of our requirement that
$g(R_{k})$ includes $q_{k}(S_{k};\xi)f(S_{k}).$

Step 2: Simplify the new estimating equation. We want to show that

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[W_{k}\left(Y_{k+\ell}-g(R_{k})^{\prime}\alpha-
A_{k}f(S_{k})^{\prime}\beta\right)(A_{k}-q_{k}(S_{k};\xi))f(S_{k})|R_{k},T_{k+\ell}<\infty\right]$
(1)
$\displaystyle=q_{k}(S_{k};\xi)(1-q_{k}(S_{k};\xi))f(S_{k})\left(\zeta(k,\ell,R_{k})-f(S_{k})^{\prime}\beta\right)+\epsilon_{k}$

for some (error) term $\epsilon_{k}$. To that end, let $Z_{k}$ be shorthand
for $Y_{k+\ell}-g(R_{k})^{\prime}\alpha-A_{k}f(R_{k})^{\prime}\beta$, and
$\rho_{k}(R_{k})$ shorthand for $\mathop{\rm
I\kern-1.99997ptE}\nolimits[A_{k}|R_{k},T_{k+\ell}<\infty].$ Notice that

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[W_{k}Z_{k}(A_{k}-q_{k}(S_{k};\xi))f(S_{k})|R_{k},T_{k+\ell}<\infty\right]$
$\displaystyle\qquad=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[W_{k}Z_{k}(A_{k}-q_{k}(S_{k};\xi))f(S_{k})|R_{k},A_{k}=1,T_{k+\ell}<\infty\right]\rho_{k}(R_{k})+$
$\displaystyle\qquad\quad\,\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[W_{k}Z_{k}(A_{k}-q_{k}(S_{k};\xi))f(S_{k})|R_{k},A_{k}=0,T_{k+\ell}<\infty\right](1-\rho_{k}(R_{k}))$
$\displaystyle\qquad=\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}|R_{k},A_{k}=1,T_{k+\ell}<\infty\right]\frac{q_{k}(S_{k};\xi)}{p_{k}(R_{k};\eta)}(1-q_{k}(S_{k};\xi))f(S_{k})\rho_{k}(R_{k})+$
$\displaystyle\quad\qquad\,\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}|R_{k},A_{k}=0,T_{k+\ell}<\infty\right]\frac{1-q_{k}(S_{k};\xi)}{1-p_{k}(R_{k};\eta)}(0-q_{k}(S_{k};\xi))f(S_{k})(1-\rho_{k}(R_{k})),$

which we re-arrange to:

$\displaystyle q_{k}(S_{k};\xi)(1-q_{k}(S_{k};\xi))f(S_{k})\left(\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}|R_{k},A_{k}=1,T_{k+\ell}<\infty\right]-\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}|R_{k},A_{k}=0,T_{k+\ell}<\infty\right]\right)+\epsilon_{k}$

where $\epsilon_{k}$ is the remaining term:

$\displaystyle
q_{k}(S_{k};\xi)(1-q_{k}(S_{k};\xi))f(S_{k})(\rho_{k}(R_{k})-p_{k}(R_{k};\eta))\left(\frac{\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}|R_{k},A_{k}=1,T_{k+\ell}<\infty\right]}{p_{k}(R_{k};\eta)}+\frac{\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}|R_{k},A_{k}=0,T_{k+\ell}<\infty\right]}{1-p_{k}(R_{k};\eta)}\right).$

Focusing on

$\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}\big{|}R_{k},A_{k}=1,T_{k+\ell}<\infty\right]-\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Z_{k}\big{|}R_{k},A_{k}=0,T_{k+\ell}<\infty\right],$

we use the definition of $Z_{k}$ and Lemma 1 to simplify this expression to

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}\big{|}R_{k},A_{k}=1,T_{k+\ell}<\infty\right]-\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[Y_{k+\ell}\big{|}R_{k},A_{k}=0,T_{k+\ell}<\infty\right]-f(R_{k})^{\prime}\beta$
$\displaystyle=\zeta(k,\ell,R_{k})-f(R_{k})^{\prime}\beta.$

Upon substituting this expression into the one above, we arrive at Equation 1.

Step 3: Show $\epsilon_{k}=0$. We want to show that $\epsilon_{k}$ is zero
under the conditions of the theorem. We consider the two cases. First assume
the conditional probability model in the denominator of weights (Equation 6)
is correct for $\theta_{\infty}\in\Theta$ so that
$\rho_{k}(R_{k})=p_{k}(R_{k};\eta_{\infty}).$ Then, it is clear that
$\epsilon_{k}$ evaluated at $\theta_{\infty}$ is zero. Alternatively, assume
the outcome model is correct (Equation 4) for $\theta_{\infty}\in\Theta$ so
that

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}|R_{k},A_{k},T_{k+\ell}<\infty]=g(R_{k})^{\prime}\alpha_{\infty}+A_{k}f(S_{k})^{\prime}\beta_{\infty}.$

Then,

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits[Y_{k+\ell}-g(R_{k})^{\prime}\alpha_{\infty}-A_{k}f(S_{k})^{\prime}\beta_{\infty}|R_{k},A_{k},T_{k+\ell}<\infty]=0$

and again, it is clear that $\epsilon_{k}$ evaluated at $\theta_{\infty}$ is
zero.

Step 4: Show $\mathbb{G}\beta_{\infty}=g$. From above, we know the conditions
of the theorem imply that, if we were to evaluate

$\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[W_{k}\left(Y_{k+\ell}-g(R_{k})^{\prime}\alpha-
A_{k}f(R_{k})^{\prime}\beta\right)(A_{k}-q_{k}(S_{k};\xi))f(S_{k})\,\big{|}\,T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty)=0$

at $\theta_{\infty}$ and apply the law of iterated expectation conditioning on
$R_{t}$, we would get

$\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[q_{k}(S_{k};\xi_{\infty})(1-q_{k}(S_{k};\xi_{\infty}))\left(\zeta(k,\ell,R_{k})-f(S_{k})^{\prime}\beta_{\infty}\right)f(S_{k})\big{|}\,T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty)=0.$

Hence,

$\displaystyle\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[q_{k}(S_{k};\xi_{\infty})(1-q_{k}(S_{k};\xi_{\infty}))f(S_{k})\zeta(k,\ell,R_{k})\big{|}\,T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty)$
$\displaystyle=\left(\sum_{k=1}^{\infty}\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[q_{k}(S_{k};\xi_{\infty})(1-q_{k}(S_{k};\xi_{\infty}))f(S_{k})f(S_{k})^{\prime}\big{|}\,T_{k+\ell}<\infty\right]\mathop{\rm
I\kern-1.99997ptP}\nolimits(T_{k+\ell}<\infty)\right)\beta_{\infty}.$

This last term is $\mathbb{G}\beta_{\infty}$ where $\mathbb{G}$ is given in
the statement of the theorem. One more application of law of iterated
expectation, this time conditioning on $S_{t}$, and the definition of
$\zeta_{\rm{marg}}(k,\ell,S_{k})$ means

$\displaystyle\mathop{\rm
I\kern-1.99997ptE}\nolimits\left[q_{k}(S_{k};\xi_{\infty})(1-q_{k}(S_{k};\xi_{\infty}))f(S_{k})\zeta(k,\ell,R_{k})\Big{|}\,T_{k+\ell}<\infty\right]$
$\displaystyle=\mathop{\rm