Incorporating nonparametric methods for estimating causal excursion effects in mobile health with zero-inflated count outcomes

ABSTRACT

In mobile health, tailoring interventions for real-time delivery is of paramount importance. Micro-randomized trials have emerged as the “gold-standard” methodology for developing such interventions. Analyzing data from these trials provides insights into the efficacy of interventions and the potential moderation by specific covariates. The “causal excursion effect,” a novel class of causal estimand, addresses these inquiries. Yet, existing research mainly focuses on continuous or binary data, leaving count data largely unexplored. The current work is motivated by the Drink Less micro-randomized trial from the UK, which focuses on a zero-inflated proximal outcome, i.e., the number of screen views in the subsequent hour following the intervention decision point. To be specific, we revisit the concept of causal excursion effect, specifically for zero-inflated count outcomes, and introduce novel estimation approaches that incorporate nonparametric techniques. Bidirectional asymptotics are established for the proposed estimators. Simulation studies are conducted to evaluate the performance of the proposed methods. As an illustration, we also implement these methods to the Drink Less trial data.

causal excursion effect, count outcome, micro-randomized trial, mobile health, structural nested mean model

1 INTRODUCTION

Mobile health (mHealth) interventions, particularly those utilizing text messages or push notifications, have been developed to promote health-related behaviors. These interventions exhibit significant potential across a wide range of health issues, from alcohol consumption reduction (Song et al., 2019) to physical activity maintenance (Lee et al., 2019). Advancements in mobile and sensing technologies now facilitate real-time tracking of an individual’s internal state and context, offering timely and personalized support (Kumar et al., 2013; Spruijt-Metz and Nilsen, 2014). This has given rise to the concept of the just-in-time adaptive intervention (JITAI), which tailors treatment in response to the evolving needs and situations of the individual, bearing the goal of delivering the right treatment on the right occasion (Nahum-Shani et al., 2018). For instance, evenings are recognized as a high-risk window for individuals with a history of excessive alcohol consumption (Day et al., 2014; Bell et al., 2020). Moreover, the burden on users due to overly frequent notifications can result in disengagement (Bell et al., 2020). Therefore, JITAIs for alcohol use could leverage evenings and the history of notification delivery as tailoring variables to determine the likelihood of delivering an intervention at the current moment.

The micro-randomized trial (MRT) has emerged as the touchstone methodology for devising these interventions (Klasnja et al., 2015; Liao et al., 2016; Qian et al., 2022; Liu et al., 2023). Within an MRT, participants undergo sequential randomization, aligning them with one of the intervention options across hundreds or even thousands of decision points. Typically, data analysis from MRTs seeks to answer 3 critical scientific questions: (1) which interventions can impact the proximal outcome; (2) in which time-varying context should the intervention be delivered; and (3) does the treatment effect change with time? The “causal excursion effect,” a novel class of causal estimand, provides solutions to these inquiries (Boruvka et al., 2018; Qian et al., 2021a). Notably, this effect can be perceived as a marginal generalization of the treatment “blip” in the structural nested mean model (Robins, 1989; 1994), as it is conditional on a few selected variables instead of all past observed variables.

Standard methods like generalized estimating equations (GEE) (Liang and Zeger, 1986) or random-effects models (Laird and Ware, 1982) may not provide consistent treatment effect estimates when data includes time-varying treatments and confounders (Sullivan Pepe and Anderson, 1994; Robins and Hernan, 2008). Structural nested mean models (Robins, 1994) and G-estimation address this by modeling causal effects of time-varying treatments on time-varying outcomes. Boruvka et al. (2018) and Qian et al. (2021a) innovatively developed weighted and centered least square estimators for estimating causal excursion effects, but focused on continuous or binary data, leaving count data unexplored. Our study is motivated by the Drink Less trial (Bell et al., 2020) where the proximal outcome, the number of screen views following notifications, is zero-inflated count data, invalidating common distributional assumptions.

Building upon Yu et al. (2023), we propose new methods for estimating causal excursion effects under potentially zero-inflated count outcomes, incorporating nonparametric nuisance estimation. Addressing technical errors in collecting randomization probabilities (Shi and Dempsey, 2023), we propose doubly robust estimating equations applicable to both MRTs and observational mHealth studies with non-randomized treatments. Our contributions are 2-fold: (1) advocating nonparametric techniques like generalized additive models and two-part models (Yu et al., 2023) for nuisance estimation under zero-inflation; and (2) proposing doubly robust estimators with bidirectional asymptotics (Yu et al., 2023) as the sample size or the number of decision points goes to infinity.

The rest of this article is organized as follows. Section 2 provides an overview of the Drink Less trial as a motivating example. The notation and causal estimands are reviewed in Section 3. Section 4 details the incorporation of nonparametric methods for estimating the causal excursion effect. Section 5 provides the data generation procedure, settings, and results of the simulation study. Next, we apply the proposed method to analyzing the Drink Less data in Section 6. We conclude with a discussion in Section 7.

2 MOTIVATING EXAMPLE: DRINK LESS

Drink Less is a behavior change app that aims to help the general adult population in the UK who want to reduce hazardous and harmful alcohol consumption (Garnett et al., 2019; 2021; Bell et al., 2020). A 30-day MRT with 349 participants was conducted to improve the push notification strategy in Drink Less (Bell et al., 2020; 2023).

Participants with a baseline Alcohol Use Disorders Identification Test (AUDIT) score of 8 or higher (Saunders et al., 1993; Bohn et al., 1995), who resided in the UK, were at least 18 years old, and desired to drink less, were recruited into the trial. Every day at 8 pm, during the trial, participants were randomly given 1 of 3 intervention options: no notification, the standard notification, or a message randomly selected from a new message bank. For more details about the messages, see Bell et al. (2020). Each option was assigned according to a static randomization probability of 40%, 30%, and 30%, respectively. In this study, we employ the total number of screen views between 8 pm and 9 pm as the proximal outcome to measure the depth of engagement with the app. Prior research has validated the use of screen views as a metric for near-term engagement (Radin et al., 2018; Perski et al., 2019). Along with the interventions, data on age, AUDIT score, gender, and other covariates were also recorded. The main goal is to explore the effects of push notifications on user engagement, whether or not these effects differ according to the user’s context, and how these effects change over time.

We first examine the distribution of the proximal outcome from the Drink Less trial. Figure 1 displays the histogram and Q–Q plot of the proximal outcome, which is the number of screen views between 8 pm and 9 pm following no notification, a new notification, or the standard notification. The results suggest that the proximal outcome is not normally distributed and is highly zero-inflated.

FIGURE 1

Graphical inspection of the distribution and normality check of the number of screen views from 8 pm to 9 pm after each type of notification. A: Histogram; B: Q-Q plot.

3 CAUSAL EXCURSION EFFECT: A REVIEW

3.1 Notations

Consider a setting with longitudinal data spanning T decision points for n participants. For each participant, the treatment assignment at time t is represented by |$A_t$|⁠. We simplify by considering only binary treatments |$A_t \in \lbrace 0,1\rbrace$|⁠, where 1 indicates the administration of treatment, and 0 indicates its absence. Detailed discussions on extensions to scenarios with multi-category treatments can be found in Web Appendix A. Let |$X_t$| denote individual and contextual information or covariates collected after time |$t-1$| and up to time t. This includes prior treatments, proximal outcomes, and the participant’s availability status |$I_t \in \lbrace 0,1\rbrace$|⁠. The value 1 implies availability for treatment at t, and 0 denotes the opposite.

After providing treatment at time t, we observe a proximal outcome |$Y_{t,\Delta }$|⁠. This outcome is a deterministic function of data collected over an interval of length |$\Delta$|⁠. For example, in the Drink Less study (Section 2), the decision time is daily at 8 pm, with the proximal outcome being screen views between 8 pm and 9 pm, resulting in |$\Delta = 1$|⁠. Other mHealth studies might have larger |$\Delta$| values with each minute as a decision point (Battalio et al., 2021). We focus on situations where |$\Delta = 1$| and assume |$Y_{t,1}$| is zero-inflated count data.

In our notation, an overbar signifies a sequence of random variables; for instance, |$\overline{A}_t$| encompasses the series |$(A_1,\cdots ,A_t)$|⁠. The data collected until time t is represented by the history |$H_t = (\overline{X}_t, \overline{A}_{t-1})$|⁠.

3.2 Causal excursion effect

In the following, we introduce the potential outcomes framework (Rubin, 1974; Robins, 1989) to define the causal excursion effect. Let |$X_t(\overline{a}_{t-1})$| be the potential covariates that would have been collected, and |$A_t(\overline{a}_{t-1})$| the treatment that would have been assigned, had the participant received the treatment sequence |$\overline{a}_{t-1} \in \lbrace 0,1\rbrace ^{t-1}$|⁠. Additionally, denote by |$Y_{t,1}(\overline{a}_{t})$| the potential proximal outcome that would have been observed had that participant received the treatment sequence |$\overline{a}_{t} \in \lbrace 0,1\rbrace ^{t}$|⁠. Here, treatments |$A_t$| and covariates |$X_t$| are also expressed as potential outcomes of previous treatment to mimic mHealth settings where covariates and treatment assignments can depend on previous treatments. The potential history at time t is represented by |$H_t(\overline{A}_{t-1})=\lbrace X_1,A_1,X_2(A_1),A_2,X_3(\overline{A}_2),\cdots ,X_t(\overline{A}_{t-1})\rbrace$|⁠.

As in Boruvka et al. (2018) and Qian et al. (2021a), we are interested in estimating the causal excursion effect of treatment |$A_t$| on |$Y_{t,1}$|⁠:

$$\begin{eqnarray} \log \frac{\mathbb {E}\left\lbrace Y_{t, 1}\left(\overline{A}_{t-1}, 1\right) \mid S_t\left(\overline{A}_{t-1}\right) = s_t, I_t\left(\overline{A}_{t-1}\right)=1\right\rbrace }{\mathbb {E}\left\lbrace Y_{t, 1}\left(\overline{A}_{t-1}, 0\right) \mid S_t\left(\overline{A}_{t-1}\right) = s_t, I_t\left(\overline{A}_{t-1}\right)=1\right\rbrace }. \end{eqnarray}$$

(1)

Here, |$S_t (\overline{A}_{t-1}) \in \mathbb {R}^p$| denotes a vector of potential moderators formed from |$H_t(\overline{A}_{t-1})$|⁠. To accommodate possibly zero-inflated count outcomes, we opt for the logarithm of the ratio of, rather than the difference in, the expected outcomes. The effect contrasts 2 excursions from the treatment protocol before time t and characterizes the treatment effect in the short term. Specifically, it assesses the effect of following the protocol until time |$t-1$| and then deviating at time t to assign treatment 1, compared to a deviation that assigns treatment 0 at the same point. By conditioning on |$I_t(\overline{A}_{t-1})=1$| and |$S_t(\overline{A}_{t-1})$|⁠, the effect is defined for only individuals available for treatment at time t and is marginalized over variables in |$H_t(\overline{A}_{t-1})$| that are not in |$S_t(\overline{A}_{t-1})$|⁠. It is critical to recognize that the causal excursion effect depends on the treatment assignment protocol |$\overline{A}_{t-1}$|⁠, diverging from traditional causal inference literature (Robins, 1994). This feature is particularly valuable in the context of real-life trials, where the impact of past treatments on user burden or disengagement can lead to significantly varied treatment effects at time t (Qian et al., 2021b). Consequently, the causal excursion effect provides a framework for improving the current treatment protocol. Setting |$S_t(\overline{A}_{t-1})$| to an empty set allows for the assessment of an intervention’s direct effect on proximal outcomes. Incorporating moderators such as days since download, treatment history, or previous outcomes further enables the evaluation of treatment effect moderation. These insights are pivotal for future mobile health application development and treatment policy formulation (Luckett et al., 2020).

When choosing |$S_t(\overline{A}_{t-1}) = H_t(\overline{A}_{t-1})$|⁠, we get the fully conditional version of the causal excursion effect, which can be expressed as

$$\begin{eqnarray} \log \frac{\mathbb {E}\left\lbrace Y_{t, 1}\left(\overline{A}_{t-1}, 1\right) \mid H_t\left(\overline{A}_{t-1}\right), I_t\left(\overline{A}_{t-1}\right)=1\right\rbrace }{\mathbb {E}\left\lbrace Y_{t, 1}\left(\overline{A}_{t-1}, 0\right) \mid H_t\left(\overline{A}_{t-1}\right), I_t\left(\overline{A}_{t-1}\right)=1\right\rbrace }. \end{eqnarray}$$

(2)

The fully conditional effect closely parallels the treatment blips in the structural nested mean model (Robins, 1994). However, our focus is solely on the immediate effect of a time-varying treatment, not the cumulative effect of all previous treatments.

3.3 Identification

To estimate the causal excursion effect from the observed data, we state 3 fundamental assumptions in causal inference.

 

Assumption 1 (Consistency)

The observed data is equal to the potential outcome under the observed treatment sequence, i.e., for |$2 \le t \le T$|⁠, |$X_t = X_t(\overline{A}_{t-1})$|⁠, |$A_t=A_t(\overline{A}_{t-1})$|⁠, |$Y_{t,1}=Y_{t,1}(\overline{A}_{t})$|⁠.

 

Assumption 2 (Positivity)

If the joint density |$\operatorname{Pr}(H_t=h_t,I_t=1)\gt 0$|⁠, then |$\operatorname{Pr}(A_t=a_t|H_t=h_t,I_t=1)\gt 0$| almost everywhere for |$a_t \in \lbrace 0,1\rbrace$|⁠.

 

Assumption 3 (Sequential ignorability)

For |$1\le t \le T$|⁠, the potential outcomes |$\lbrace X_{t+1}(\overline{a}_t),A_{t+1}(\overline{a}_t), \ldots , X_{T+ 1}(\overline{a}_T)\rbrace$| are independent of |$A_t$| conditional on |$H_t$|⁠.

Assumption 1 connects the potential outcomes with the observed data and states that there is no interference between the observations. Both Assumptions 2 and 3 are inherently satisfied in MRTs due to the sequential randomization of treatments based on known probabilities.

Under these assumptions, the causal excursion effect can be expressed in terms of observed data:

$$\begin{eqnarray} \log \frac{\mathbb {E}\left[\mathbb {E}\left\lbrace Y_{t, 1} \mid A_t=1, H_t, I_t=1\right\rbrace \mid S_t, I_t=1\right]}{\mathbb {E}\left[\mathbb {E}\left\lbrace Y_{t, 1} \mid A_t=0, H_t, I_t=1\right\rbrace \mid S_t, I_t=1\right]}. \end{eqnarray}$$

(3)

The derivation of the identifiability results follows similarly to Qian et al. (2021a). For completeness, the proof for the general case where |$\Delta \gt 1$| is included in Web Appendix B.

Additionally, the causal excursion effect simplifies to

$$\begin{eqnarray} \log \frac{\mathbb {E}\left(Y_{t, 1} \mid A_t=1, H_t, I_t=1\right)}{\mathbb {E}\left(Y_{t, 1} \mid A_t=0, H_t, I_t=1\right)} . \end{eqnarray}$$

(4)

When conditioning on the full history |$H_t$|⁠. In what follows, we begin by detailing the estimation procedure for the fully conditional excursion effect, subsequently extending the approach to instances where |$S_t\ne H_t$|⁠.

4 INCORPORATING NONPARAMETRIC METHODS FOR CAUSAL EXCURSION EFFECT ESTIMATION

4.1 Estimating the conditional excursion effect

In this section, we discuss the estimation procedure for the fully conditional excursion effect, drawing inspiration from the G-estimator initially introduced by Robins (1994) and subsequently discussed by Yu et al. (2023) for the multiplicative structural nested mean model. Our approach differs in several aspects: it centers on estimating the causal excursion effect, which serves as a marginal extension of the treatment “blip” in structural nested mean models, as detailed in Section 4.2. Additionally, it introduces a weighting component to enable semi-parametric efficiency or facilitate the estimation of marginal treatment effects. The addition of an availability indicator, |$I_t$|⁠, further adapts our model for effective application in mHealth contexts.

We begin by assuming a parametric model for the conditional excursion effect:

$$\begin{eqnarray} \log \frac{\mathbb {E}\left(Y_{t, 1} \mid A_t=1, H_t, I_t=1\right)}{\mathbb {E}\left(Y_{t, 1} \mid A_t=0, H_t, I_t=1\right)} = f(H_t)^{\top }\phi \end{eqnarray}$$

(5)

for unknown p-dimensional parameter vector |$\phi$|⁠, where |$f(\cdot )$| is a known deterministic function. Throughout, we denote the true value of |$\phi$| by |$\phi _0$|⁠.

Then, define |$U_t(\phi ) = Y_{t,1} \exp (-A_tf(H_t)^{\top }\phi )$| to mimic the potential outcome |$Y_{t,1}(\overline{A}_{t-1},0)$| that would have been observed had the treatment been removed at time t. By Assumption 3,

$$\begin{eqnarray} \mathbb {E}[U_t(\phi )|A_t, H_t, I_t=1] &=& \mathbb {E}[Y_{t,1}(\overline{A}_{t-1},0)|A_t,H_t, I_t=1] \\ &=& \mathbb {E}[Y_{t,1}(\overline{A}_{t-1},0)|H_t, I_t=1] \\ &=& \mathbb {E}[U_t(\phi )|H_t, I_t=1]. \end{eqnarray}$$

(6)

Let |$h_0(H_t) = \mathbb {E}[U_t(\phi )|H_t, I_t=1]$| denote the conditional mean of |$U_t(\phi )$| given |$H_t$| and |$I_t=1$|⁠, and let |$p_{t,0}(H_t)=\operatorname{Pr}(A_t=1|H_t)$| denote the treatment randomization probability. Estimation of |$\phi$| in model (5) can thus be based on the following estimating function:

$$\begin{eqnarray} m_C(H_{t};\phi , h, p_t) &=& I_t \left[U_t(\phi ) - h(H_t)\right]\\ &&\times \, \left[\lbrace A_t - p_t(H_t)\rbrace f(H_t)\right], \end{eqnarray}$$

(7)

for any given h and |$p_t$|⁠. Building upon argument (6), the expectation of (7) equals zero when |$\phi =\phi _0$| under conditions |$h = h_0$| or |$p_t = p_{t,0}$|⁠. To denote individual participants, we introduce the subscript i, resulting in n independent and identically distributed copies of the data sequence |$H_T$|⁠, represented as |$H_{1,T}, \ldots , H_{n,T}$|⁠. The estimating equation for |$\phi$| is expressed as

$$\begin{eqnarray} \frac{1}{nT}\sum _{i=1}^n\sum _{t=1}^T m_C(H_{it};\phi , h, p_t) = 0. \end{eqnarray}$$

(8)

In MRTs with known randomization probabilities, estimating the conditional mean |$h_0(H_t)$| is the primary task. We express it as

$$\begin{eqnarray} h_0(H_t) &=& \mu _{1,0}(H_t)\exp \lbrace -f(H_t)^{\top }\phi \rbrace p_{t}(H_t) \\ &&+\, \mu _{0,0}(H_t)\left(1-p_{t}(H_t)\right), \end{eqnarray}$$

where |$\mu _{1,0}(H_t) = \mathbb {E}[Y_{t,1}|A_t=1,H_t, I_t=1]$| and |$\mu _{0,0}(H_t) = \mathbb {E}[Y_{t,1}|A_t=0,H_t, I_t=1]$|⁠. Nonparametric models such as generalized additive models can be employed to estimate |$\mu _{1,0}(H_t)$| and |$\mu _{0,0}(H_t)$|⁠. Plugging estimates |$\widehat{\mu }_1$| and |$\widehat{\mu }_0$| into (7) yields the estimating function |$m_C(H_t; \phi , \widehat{\mu }_1, \widehat{\mu }_0, p_{t,0})$|⁠. Denote |$\eta =(\mu _0, \mu _1, p_t)$|⁠, the true value |$\eta _0 = (\mu {0,0}, \mu {1,0}, p{t,0})$|⁠, and its estimator |$\widehat{\eta }=(\widehat{\mu }_0, \widehat{\mu }_1, \widehat{p}_t)$|⁠. For simplicity, denote |$m_C(H_t; \phi , \eta ) = m_C(H_t; \phi , \mu _1, \mu _0, p_t)$|⁠. This resulting estimator is referenced as estimator for the conditional effect that incorporates nonparametric methods (ECE-NonP), following the terminology used by Qian et al. (2021a).

In observational studies, both |$h_0(H_t)$| and |$p_{t,0}(H_t)$| are unknown, necessitating estimation of these nuisance functions from data. Since we cannot guarantee consistent nuisance estimators, a doubly robust estimator ensuring consistency when either |$h_0(H_t)$| or |$p_{t,0}(H_t)$| is consistently estimated becomes crucial. We demonstrate the rate double robustness of our proposed ECE-NonP estimator below.

Following Yu et al. (2023), we make an additional assumption:

 

Assumption 4

As |$n T \rightarrow \infty ,\Vert \widehat{\eta }-\eta _0\Vert _{2, P} \rightarrow 0$| in probability and

$$\begin{eqnarray} &&\left(\left\Vert \widehat{\mu }_1-\mu _{1,0}\right\Vert _{2, P}+\left\Vert \widehat{\mu }_0-\mu _{0,0}\right\Vert _{2, P}+\left\Vert \widehat{p}_t-p_{t,0}\right\Vert _{2, P}\right)\\ &&\quad\Vert \widehat{p}_t- p_{t,0}\Vert _{2, P}=o_P\left\lbrace (n T)^{-1 / 2}\right\rbrace , \end{eqnarray}$$

which outlines nuisance function convergence rates for establishing asymptotics. The first part ensures estimator consistency, while the second part determines the ECE-NonP asymptotic distribution. For MRT data, this assumption is naturally met as |$p_{t,0}(H_t)$| is known, and others can be consistently estimated using nonparametric methods. Theorem 1 provides bidirectional asymptotics for ECE-NonP, with the proof presented in Web Appendix C.

 

Theorem 1 (Bidirectional Asymptotics of ECE-NonP)

Suppose model (5), Assumptions 1–3, Assumption 4, and some regularity conditions hold. Let |${\dot{m}}_\phi (H_t; \phi , \eta ) = \frac{\partial m_C(H_t; \phi ,\eta )}{\partial \phi }$|⁠. Let |$\mathbb {P}g\lbrace m_C(H_{t}; \phi , \eta )\rbrace = \frac{1}{T}\sum _{t=1}^{T} \mathbb {E}[g\lbrace m_C(H_{i,t};\phi , \eta )\rbrace ]$|⁠, where g is any given function or operator of |$m_{C}$|⁠. As either |$n \rightarrow \infty$| or |$T \rightarrow \infty$|⁠, we have |$(nT)^{1/2}(\widehat{\phi} - \phi_0) \rightarrow \text{MVN}(0, \boldsymbol {B} \boldsymbol {\Sigma } \boldsymbol {B}^{\top })$| in distribution, where |$\boldsymbol {B} = \lbrace \mathbb {P}{\dot{m}}_\phi (H_{t};\phi _0,\eta _0)\rbrace ^{-1}$| and

$$\begin{eqnarray} \boldsymbol{\Sigma} &=& \lim_{n T \rightarrow \infty} \frac{1}{nT} \sum _{i=1}^{n} \sum _{t=1}^{T}\\ && \mathbb {E}\left\lbrace m_{C}(H_{i,t};\phi _0, \eta _0) m_{C}^{\top}(H_{i,t}; \phi _0, \eta _0)|H_{i,t}\right\rbrace . \end{eqnarray}$$

 

Remark 1

Qian et al. (2021a) introduced an estimator ECE for the conditional excursion effect with binary outcomes. Following Kim et al. (2021), ECE can be extended to count outcomes using the estimating function:

$$\begin{eqnarray} &&I_t \exp \lbrace -A_t f(H_t)^{\top } \phi \rbrace [Y_{t, 1}-\exp \lbrace g(H_t)^{\top } \alpha +A_t f(H_t)^{\top } \phi \rbrace]\\ &&\quad \times \widetilde{K}_t\begin{bmatrix}g(H_t) \\ \left\lbrace A_t-p_t(H_t)\right\rbrace f(H_t) \end{bmatrix}, \end{eqnarray}$$

where |$\widetilde{K}_t$| is defined as

$$\begin{eqnarray} -\frac{\exp \left\lbrace f(H_t)^{\mathrm{T}} \phi \right\rbrace }{\exp \left\lbrace f(H_t)^{\mathrm{T}} \phi \right\rbrace p_t(H_t)+\left\lbrace 1-p_t(H_t)\right\rbrace }. \end{eqnarray}$$

Their estimator assumes a parametric working model for |$\mathbb {E}(Y_{t,1}|H_t, A_t=0, I_t=1)$|⁠. Instead, we represent |$h_0(H_t)$| using nonparametrically estimated |$\widehat{\mu }_{1}(H_t)$| and |$\widehat{\mu }_{0}(H_t)$|⁠. Semiparametric efficiency is achieved when the working model is correctly specified. We examine ECE and ECE-NonP’s performance with the above weight |$\widetilde{K}_t$| in simulations.

 

Remark 2

The ECE-NonP estimator’s consistency relies on correctly specifying the treatment effect model (5). However, the high dimensionality of the full history |$H_t$| complicates accurate specification of the model. Currently, there is limited guidance on how to accurately determine the treatment effect model, which remains a prominent challenge.

 

Remark 3

The principle of bidirectional asymptotics suggests that as either the sample size n or the number of decision points T approaches infinity, the ECE-NonP estimator achieves consistency and asymptotic normality. However, this should be used with caution. A larger T may not necessarily imply more information for estimating |$\phi$| when considering user availability during a trial. Time points where |$I_t=0$| do not contribute to (7), and the value of |$I_t$| can vary over time. This is not a problem in studies like Drink Less, where users are assumed to be always available given the intervention nature.

4.2 Estimating the marginal excursion effect

For the marginal excursion effect, we assume the following parametric model for |$1 \le t \le T$|⁠:

$$\begin{eqnarray} \log \frac{\mathbb {E}\left[\mathbb {E}\left\lbrace Y_{t, 1} \mid A_t=1, H_t, I_t=1\right\rbrace \mid S_t, I_t=1\right]}{\mathbb {E}\left[\mathbb {E}\left\lbrace Y_{t, 1} \mid A_t=0, H_t, I_t=1\right\rbrace \mid S_t, I_t=1\right]} = S_t^{\top }\beta, \\ \end{eqnarray}$$

(9)

where |$\beta$| is the parameter of interest with true value |$\beta _0$|⁠. A linear model is assumed for ease of interpretation, but the procedure extends to nonlinear forms (Qian et al., 2021a).

Then, we construct a variable |$U_t = Y_{t,1}\exp (-A_tS_t^{\top }\beta )$| under model (9). Estimation of |$\beta$| can be based on the following estimating function:

$$\begin{eqnarray} m_M(H_t; \beta , \eta ) = I_t W_t\left[U_t(\beta ) - h(H_t)\right]\times \left[\lbrace A_t - \widetilde{p}_t(S_t)\rbrace S_t\right], \\ \end{eqnarray}$$

(10)

where weight |$W_t$| is added to allow for the estimation of marginal causal excursion effects conditional on |$S_t$| instead of |$H_t$|⁠. The weight |$W_t$| is defined as

$$\begin{eqnarray} \left\lbrace \frac{\widetilde{p}_t(S_t)}{p_t(H_t)}\right\rbrace ^{A_t}\left\lbrace \frac{1-\widetilde{p}_t(S_t)}{1-p_t(H_t)}\right\rbrace ^{1-A_t} , \end{eqnarray}$$

where the numerical probability |$\widetilde{p}_t(S_t)$| can be chosen arbitrarily as long as it only depends on |$S_t$|⁠. Intuitively, the weight |$W_t$| transforms the data distribution where |$A_t$| is randomized with probability |$p_t(H_t)$| to a distribution where |$A_t$| is randomized with probability |$\widetilde{p}_t(S_t)$|⁠. The solution to |$\frac{1}{nT}\sum _{i=1}^n\sum _{t=1}^{T}m_M(H_{i,t}; \beta , \eta ) = 0$| yields an estimator for |$\beta$|⁠.

Similar to before, we need to estimate the nuisance function |$h_0(H_t)$| and plug estimates into (10). We can express |$h_0(H_t)$| using 2 conditional means:

$$\begin{eqnarray} h_0(H_t) &=& \mu _{1,0}(H_t) \exp \lbrace -S_t^{\top }\beta \rbrace \widetilde{p}_t(S_t)\\ && +\, \mu _{0,0}(H_t)\left(1-\widetilde{p}_t(S_t)\right), \end{eqnarray}$$

where |$p_t(H_t)$| is replaced by |$\widetilde{p}_t(S_t)$| due to the change of probability induced by weight |$W_t$|⁠. In MRTs, we assume the propensity score |$p_t(H_t)$| appearing in |$W_t$| is known. This estimator is referred to as the estimator for the marginal excursion effect that incorporates nonparametric methods (EMEE-NonP).

Theorem 2 provides the bidirectional asymptotics of the EMEE-NonP estimator, proven in Web Appendix D.

 

Theorem 2 (Bidirectional Asymptotics of EMEE-NonP)

Suppose the randomization probability |$p_t(H_t)$| is known. Suppose model (9), Assumptions 1–3, Assumption 4, and some regularity conditions hold. Let |${\dot{m}}_\beta (H_t; \beta , \eta ) = \frac{\partial m_M(H_t; \beta ,\eta )}{\partial \beta }$|⁠. Let |$\mathbb {P}g\lbrace m_M(H_{t}; \beta , \eta )\rbrace = \frac{1}{T}\sum _{t=1}^{T}\mathbb {E}[g\lbrace m_M(H_{i,t};\beta , \eta )\rbrace ]$| where g is any given function or operator of |$m_{M}$|⁠. As either |$n \rightarrow \infty$| or |$T \rightarrow \infty$|⁠, we have |$(nT)^{1/2}(\widehat{\beta } - \beta _0) \rightarrow \operatorname{MVN}(0, \boldsymbol {B} \boldsymbol {\Sigma } \boldsymbol {B}^{\top })$| in distribution, where |$\boldsymbol {B} = \lbrace \mathbb {P}{\dot{m}}_\beta (H_{t};\beta _0,\eta _0)\rbrace ^{-1}$| and

$$\begin{eqnarray} \boldsymbol {\Sigma } &=& \lim _{n T \rightarrow \infty } \frac{1}{n T} \sum _{i=1}^n \sum _{t=1}^{T} \mathbb {E}\left\lbrace m_{M}\left(H_{i,t};\beta _0, \eta _0\right)\right.\\ &&\left. m_{M}^{\top }\left(H_{i,t}; \beta _0, \eta _0\right) \mid H_{i,t}\right\rbrace . \end{eqnarray}$$

 

Remark 4

The EMEE method, as proposed in Qian et al. (2021a) for a binary outcome, can concurrently determine both |$\alpha$| and |$\beta$|⁠. The estimating equation is expressed as

$$\begin{eqnarray} &&\sum _{t=1}^{T} I_t \exp \left\lbrace -A_t S_t^{\top } \beta \right\rbrace \left[Y_{t, 1}-\exp \left\lbrace g(H_t)^{\top } \alpha +A_t S_t^{\top } \beta \right\rbrace \right]\\ &&\times \ W_t\begin{bmatrix}g(H_t) \\ \left\lbrace A_t-\widetilde{p}_t(S_t)\right\rbrace S_t \end{bmatrix} . \end{eqnarray}$$

(11)

One important property of the estimator is that the consistency is robust to the misspecification of the working model |$g(H_t)^{\top }\alpha$|⁠. In our method, the essential requirement for robustness is embodied in Assumption 4.

 

Remark 5

The EMEE-NonP estimator lacks the double robustness property, which we demonstrate:

$$\begin{eqnarray} && \mathbb {E}\left[S_tI_tW_t\left\lbrace U_t(\beta ) - h(H_t)\right\rbrace \times \left\lbrace A_t-\widetilde{p}_t(S_t)\right\rbrace \right]\\ \quad &=& \mathbb {E}\left[\mathbb {E}\left\lbrace S_tI_tW_t\left(U_t(\beta ) {-} h(H_t)\right)\times \left(A_t-\widetilde{p}_t(S_t)\right)|A_t, H_t\right\rbrace \right]\\ \quad &=& \mathbb {E}\left\lbrace S_tI_t\mathbb {E}\left[W_t |A_t, H_t\right] \times [\mathbb {E}(U_t(\beta )|A_t, H_t, I_t=1) {-} h(H_t)]\right.\\ && \left.\times \left[p_t(H_t)-\widetilde{p}_t(S_t)\right]\right\rbrace . \end{eqnarray}$$

For double robustness, the condition |$\mathbb {E}(U_t(\beta )|A_t, H_t,I_t=1) = h_0(H_t) = \mathbb {E}[U_t(\beta )|H_t, I_t=1]$| is required. However, under Assumptions 1–  3, this cannot be ensured. A similar issue arises for the EMEE estimator. As highlighted in Qian et al. (2021b), even if |$\exp \lbrace g(H_t)^{\top } \alpha \rbrace$| correctly models |$\mathbb {E}\lbrace Y_{t, 1}(\overline{A}_{t-1}, 0) \mid H_t, I_t=1, A_t=0\rbrace$|⁠, the term |$Y_{t, 1}-\exp \lbrace g(H_t)^{\top } \alpha +A_t S_t^{\top } \beta \rbrace$| in the estimating Equation (11) does not generally have conditional expectation zero given |$H_t$|⁠. This occurs because, while the causal excursion effect exists as a marginal model, the fully conditional effect, dependent on |$H_t$|⁠, is not necessarily |$S_t^{\top } \beta$|⁠.

 

Remark 6

For the estimator of |$\beta$| to be consistent, the treatment effect model defined in (9) must be correctly specified. Since |$S_t$| consists of summary variables chosen from |$H_t$| by the researcher, this is an easier task compared to modeling the conditional excursion effect. Yet, by conditioning on only part of the history, i.e., |$S_t$|⁠, the marginal excursion effect depends on the randomization probability of past treatment assignments |$\overline{A}_{t-1}$| (Guo et al., 2021; Zhang and Laber, 2021; Qian et al., 2021b). Therefore, any interpretation of the causal excursion effect must be contextualized within the current treatment protocol.

4.3 A doubly robust EMEE-NonP estimator

The preceding EMEE-NonP estimator lacks double robustness, which means that it solely relies on known or correctly specified randomization probabilities. This makes it difficult to extend to observational mHealth studies. In this section, we provide a doubly robust estimator, referred to as the DR-EMEE-NonP estimator.

The modified estimating function is expressed as

$$\begin{eqnarray} m_D(H_t; \beta , \eta ) &=& I_t W_t\left\lbrace U_t(\beta ) - \mu _{A_t}(H_t)\exp (-A_tS_t^\top \beta )\right\rbrace \\ &&\times \, [\lbrace A_t - \widetilde{p}_t(S_t)\rbrace S_t] \\ &&+\, I_t\widetilde{p}_t(S_t)(1-\widetilde{p}_t(S_t))\\ &&\times \,\left\lbrace \mu _{1}(H_t)\exp (-S_t^{\top }\beta )-\mu _0(H_t)\right\rbrace S_t, \\ \end{eqnarray}$$

(12)

where |$\mu _{A_t}(H_t) = \mathbb {E}[Y_{t,1}|H_t,A_t,I_t=1]$|⁠.

Below we make a new assumption about the convergence rate of nuisance functions and establish the bidirectional asymptotics for DR-EMEE-NonP in Theorem 3. The proof is deferred to Web Appendix E.

 

Assumption 5

As |$n T \rightarrow \infty ,\left\Vert \widehat{\eta }-\eta _0\right\Vert _{2, P} \rightarrow 0$| in probability and

$$\begin{eqnarray} &&\left(\Vert \widehat{\mu }_1 - \mu _{1,0}\Vert _{2,P} + \Vert \widehat{\mu }_0 - \mu _{0,0}\Vert _{2,P}\right)\Vert \widehat{p}_t - p_{t,0}\Vert _{2,P} \\ &&\quad= o_P\lbrace (nT)^{-1/2}\rbrace . \end{eqnarray}$$

 

Theorem 3 (Bidirectional Asymptotics of DR-EMEE-NonP)

Suppose model (9), Assumptions 1–3, Assumption 5, and some regularity conditions hold. Let |${\dot{m}}_\beta (H_t; \beta , \eta ) = \frac{\partial m_D(H_t; \beta ,\eta )}{\partial \beta }$|⁠. Let |$\mathbb {P}g\lbrace m_D(H_{t}; \beta , \eta )\rbrace = \frac{1}{T}\sum _{t=1}^{T} \mathbb {E}[g\lbrace m_D(H_{i,t};\beta , \eta )\rbrace ]$|⁠, where g is any given function or operator of |$m_{D}$|⁠. As either |$n \rightarrow \infty$| or |$T \rightarrow \infty$|⁠, we have |$(nT)^{1/2}(\widehat{\beta } - \beta _0) \rightarrow \operatorname{MVN}(0, \boldsymbol {B} \boldsymbol {\Sigma } \boldsymbol {B}^{\top })$| in distribution, where |$\boldsymbol {B} = \lbrace \mathbb {P}{\dot{m}}_\beta (H_{t};\beta _0,\eta _0)\rbrace ^{-1}$| and

$$\begin{eqnarray} \boldsymbol {\Sigma } &=& \lim _{n T \rightarrow \infty } \frac{1}{n T} \sum _{i=1}^n \sum _{t=1}^{T}\\ &&\times \, \mathbb {E}\left\lbrace m_{D}\left(H_{i,t};\beta _0, \eta _0\right) m_{D}^{\top }\left(H_{i,t}; \beta _0, \eta _0\right) \mid H_{i,t}\right\rbrace . \end{eqnarray}$$

The idea of forming a doubly robust estimator originated with the research of Scharfstein et al. (1999). This foundational work was subsequently expanded upon by multiple studies, notably the DR-Learner as explored by van der Laan and Rubin (2006), Nie and Wager (2021), and Kennedy (2023). More recently, Shi and Dempsey (2023) introduced a DR-learner designed to estimate causal excursion effects in mHealth studies, specifically for continuous and binary outcomes. Our work distinguishes itself in 2 primary ways. Firstly, the proposed approaches center on estimating the causal excursion effect for zero-inflated count outcomes—an issue previously unexamined. Secondly, we establish bidirectional asymptotics for the DR-EMEE-NonP estimator, necessitating that either the sample size or the number of decision points approaches infinity.

5 SIMULATION STUDIES

In this section, we conduct extensive simulation experiments to assess the finite-sample performance of the proposed estimators across several scenarios:

1. A binary treatment MRT with |$p_t (H_t)=\operatorname{expit}(-0.5 A_{t-1}+ 0.5 Z_t)$|⁠, where |$Z_t$| is a time-varying covariate.

2. A binary treatment observational study with |$p_t(H_t)=\operatorname{expit}(-0.5 A_{t-1}+ 0.5 Z_t)$|⁠.

3. A binary treatment MRT with |$p_t(H_t)$| given by a Thompson sampling (TS) algorithm (Russo et al., 2018).

4. A 3-category treatment MRT with |$p_{t,1}(H_t)=p_{t,2}(H_t) =0.5\operatorname{expit}(-0.5 A_{t-1}+ 0.5 Z_t)$|⁠.

In Scenarios (1) and (2), we compare 8 estimators: ECE, ECE-NonP, EMEE, EMEE-NonP, DR-EMEE-NonP, GEE with independence (GEE.IND) or exchangeable (GEE.EXCH) working correlation, and the G-estimator (Yu et al., 2023), which mirrors ECE-NonP without |$\widetilde{K}_t$| and EMEE-NonP without |$W_t$|⁠. Other scenarios focus on EMEE, EMEE-NonP, DR-EMEE-NonP, and GEE methods.

We adopt a simple setting with |$\Delta = 1$|⁠, and all participants are available at all decision points, i.e., |$I_t=1$|⁠, |$t=1,\ldots ,T$|⁠. The time-varying covariate |$Z_t$| can take 3 values, 0, 1, and 2, each with an equal probability. We generate the outcome |$Y_{t,1}$| using a zero-inflated negative binomial (NB) model as follows:

$$\begin{eqnarray} O_{t} &\sim & \operatorname{Bernoulli}(\pi _{t}),\\ L_{t} &\sim & \operatorname{NB}(\mu _{t}, r),\\ Y_{t,1} &=& O_{t}L_{t}, \end{eqnarray}$$

where |$r=1$| denotes the dispersion parameter.

In Scenario (1), we set |$\pi _{t} = \exp (-0.4(Z_t + 0.1) + 0.1 Z_t A_t)$| and |$\mu _{t} = \lbrace 2.2 \mathbb {1}_{Z_t=0}+2.5 \mathbb {1}_{Z_t=1} +2.4 \mathbb {1}_{Z_t=2}\rbrace \exp \lbrace A_t(0.1+0.3 Z_t)\rbrace$|⁠. Hence, the true conditional causal excursion effect is

$$\begin{eqnarray} \log \frac{\mathbb {E}\left(Y_{t,1} \mid A_t=1, H_t\right)}{\mathbb {E}\left(Y_{t,1} \mid A_t=0, H_t\right)}=0.1+0.4 Z_t. \end{eqnarray}$$

Here, |$Z_t$| interacts with the treatment |$A_t$|⁠, a moderator that impacts the conditional causal excursion effect. We also consider the fully marginal excursion effect, which is

$$\begin{eqnarray} \log \frac{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t,1} \mid A_t=1, H_t\right)\right\rbrace }{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t,1} \mid A_t=0, H_t\right)\right\rbrace }=0.460. \end{eqnarray}$$

For Scenario (2), we set |$\pi _{t} = \exp (-0.4(Z_t + 0.1) + 0.1 Z_t A_t)$| and |$\mu _{t} = \exp \lbrace 0.2 + 0.5Z_t + A_t(0.1+0.3 Z_t)\rbrace$|⁠. In this case, the true marginal excursion effect becomes

$$\begin{eqnarray} \log \frac{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t,1} \mid A_t=1, H_t\right)\right\rbrace }{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t,1} \mid A_t=0, H_t\right)\right\rbrace }=0.578. \end{eqnarray}$$

We omit the details of Scenarios (3) and (4) here and defer these to Web Appendix F.

The numerator of the weight |$W_t$|⁠, or the numerical probability |$\widetilde{p}_t(S_t)$|⁠, is a constant in t, given by |$\frac{\sum _{n=1}\sum _{t=1}^TA_{i,t}}{Tn}$|⁠. A working model |$g(H_t)^{\top }\alpha = \alpha _0 + \alpha _1 Z_t$| is used for the logarithm of the expected outcome under no treatment when implementing ECE, EMEE, and GEE methods. For ECE-NonP, EMEE-NonP, and DR-EMEE-NonP, the preliminary step involves estimating the nuisance functions |$\mathbb {E}[Y_{t,1}|A_t=1,H_t]$| and |$\mathbb {E}[Y_{t,1}|A_t=0,H_t]$| via nonparametric regressions, specifically using generalized additive models. Given the zero-inflation in the data, a hurdle model is used to model the conditional mean in 2 parts: the probability of attaining value 0 and the non-zero counts. (Hu et al., 2011). In the observational data framework, the objective is to examine the rate double robustness of the proposed DR-EMEE-NonP estimator. Here, the randomization probability |$p_{t,0}(H_t)$| also requires estimation from the data, for which we use the sample proportion.

The performance measures include estimation bias (Bias), mean estimated standard error (SE), standard deviation (SD), root mean squared error (RMSE), and coverage probability of 95% confidence interval (CP) across 1000 replicates. In the simulation experiments, we set the number of decision points to |$T= 30, 100, 150$| and the sample size to |$n = 100$|⁠.

Table 1 presents simulation results of marginal excursion effects under Scenario (1). DR-EMEE-NonP, EMEE-NonP, and EMEE consistently display negligible bias across all settings, with empirical coverage probabilities of 95% CI closely aligned to the nominal level. In contrast, ECE-NonP, ECE, and the G-estimator underperform compared to DR-EMEE-NonP, EMEE-NonP, and EMEE, exhibiting greater biases and empirical coverage probabilities of 95% CI significantly deviating from the nominal level due to misspecification of the treatment effect model.

TABLE 1

Comparison of 8 estimators for the fully marginal excursion effect under Scenario (1).

Estimator	Time length	Bias	SE	SD	RMSE	CP
	30	−0.025	0.055	0.056	0.061	0.91
	100	−0.025	0.030	0.031	0.040	0.86
ECE	150	−0.024	0.025	0.025	0.035	0.83
	30	−0.025	0.056	0.056	0.061	0.92
	100	−0.025	0.030	0.031	0.040	0.87
ECE-NonP	150	−0.024	0.025	0.025	0.035	0.84
	30	−0.013	0.057	0.058	0.059	0.93
	100	−0.013	0.032	0.032	0.035	0.92
G-estimator	150	−0.012	0.026	0.026	0.029	0.92
	30	−0.001	0.058	0.059	0.059	0.94
	100	−0.001	0.032	0.033	0.033	0.94
EMEE	150	0.000	0.026	0.026	0.026	0.93
	30	−0.001	0.058	0.058	0.058	0.95
	100	−0.001	0.032	0.033	0.033	0.94
EMEE-NonP	150	0.000	0.026	0.026	0.026	0.94
	30	−0.001	0.058	0.059	0.058	0.95
	100	−0.001	0.032	0.033	0.033	0.94
DR-EMEE-NonP	150	0.000	0.026	0.026	0.026	0.94
	30	−0.026	0.055	0.055	0.061	0.91
	100	−0.025	0.030	0.031	0.040	0.87
GEE (ind)	150	−0.025	0.024	0.025	0.035	0.83
	30	−0.026	0.055	0.056	0.061	0.91
	100	−0.025	0.030	0.031	0.040	0.87
GEE (exch)	150	−0.025	0.024	0.025	0.035	0.83

Bias or coverage probability is in bold if it deviate significantly from 0 or 0.95, respectively.

Abbreviations: CP, coverage probability; GEE, generalized estimating equations; RMSE, root mean squared error; SD, standard deviation; SE, standard error.

Table 2 displays the simulation results of conditional excursion effects under Scenario (1). ECE, ECE-NonP, and the G-estimator parallel the performances of EMEE, EMEE-NonP, and DR-EMEE-NonP, exhibiting minimal bias and great empirical coverage probabilities due to the correct specification of the treatment effect model by including the moderator |$Z_t$|⁠. However, GEE methods yield biased results across all settings in Tables 1 and 2, with empirical coverage probabilities notably diverging from nominal values.

TABLE 2

Comparison of 8 estimators of the treatment effect moderation |$Z_t$| under Scenario (1).

		|$\beta _0$|					|$\beta _1$|
Estimator	Time length	Bias	SE	SD	RMSE	CP	Bias	SE	SD	RMSE	CP
	30	−0.003	0.076	0.077	0.077	0.94	0.004	0.076	0.077	0.077	0.95
	100	0.001	0.041	0.043	0.043	0.94	−0.002	0.041	0.043	0.043	0.92
ECE	150	−0.001	0.034	0.034	0.034	0.95	0.001	0.034	0.035	0.035	0.94
	30	−0.003	0.076	0.078	0.078	0.94	0.004	0.076	0.077	0.077	0.95
	100	0.001	0.042	0.043	0.043	0.95	−0.002	0.042	0.043	0.043	0.93
ECE-NonP	150	−0.001	0.034	0.034	0.034	0.95	0.001	0.034	0.035	0.035	0.95
	30	−0.003	0.076	0.077	0.077	0.94	0.004	0.076	0.077	0.077	0.95
	100	0.001	0.042	0.042	0.042	0.95	−0.002	0.041	0.043	0.043	0.93
G-estimator	150	−0.001	0.034	0.034	0.034	0.95	0.001	0.034	0.035	0.035	0.95
	30	−0.002	0.076	0.077	0.077	0.94	0.003	0.076	0.078	0.078	0.94
	100	0.001	0.041	0.042	0.042	0.94	−0.002	0.041	0.043	0.043	0.93
EMEE	150	−0.001	0.034	0.034	0.034	0.94	0.001	0.034	0.035	0.035	0.95
	30	−0.003	0.076	0.077	0.077	0.94	0.003	0.076	0.078	0.078	0.95
	100	0.001	0.042	0.042	0.042	0.95	−0.002	0.042	0.043	0.043	0.93
EMEE-NonP	150	−0.001	0.034	0.034	0.034	0.95	0.001	0.034	0.035	0.035	0.95
	30	−0.003	0.076	0.077	0.077	0.94	0.003	0.076	0.078	0.078	0.95
	100	0.001	0.042	0.042	0.042	0.95	−0.002	0.042	0.043	0.043	0.93
DR-EMEE-NonP	150	−0.001	0.034	0.034	0.034	0.95	0.001	0.034	0.035	0.035	0.95
	30	0.021	0.074	0.075	0.078	0.94	−0.024	0.072	0.073	0.077	0.92
	100	0.024	0.040	0.041	0.048	0.89	−0.027	0.039	0.041	0.049	0.88
GEE (ind)	150	0.022	0.033	0.033	0.039	0.90	−0.024	0.032	0.033	0.041	0.87
	30	0.021	0.074	0.075	0.078	0.94	−0.024	0.072	0.073	0.077	0.92
	100	0.024	0.040	0.041	0.048	0.89	−0.027	0.039	0.041	0.049	0.88
GEE (exch)	150	0.022	0.033	0.033	0.039	0.90	−0.025	0.032	0.033	0.041	0.87

Bias or coverage probability is in bold if it deviates significantly from 0 or 0.95, respectively.

Abbreviations: CP, coverage probability; GEE, generalized estimating equations; RMSE, root mean squared error; SD, standard deviation; SE, standard error.

Table 3 summarizes simulation results of marginal excursion effects under Scenario (2). Only DR-EMEE-NonP performs well, exhibiting minimal bias and impressive empirical coverage probabilities due to its rate double robustness property lacking in other methods. Yet, in Table 4, all methods exhibit excellent performances attributed to correct specification of the outcome model and treatment effect model. Scenarios (3) and (4) yield similar observations; details are provided in Web Appendix G.

TABLE 3

Comparison of 8 estimators for the fully marginal excursion effect under Scenario (2).

Estimator	Time length	Bias	SE	SD	RMSE	CP
	30	−0.015	0.065	0.067	0.068	0.93
	100	−0.015	0.036	0.036	0.039	0.92
ECE	150	−0.015	0.029	0.030	0.034	0.91
	30	−0.017	0.062	0.067	0.069	0.92
	100	−0.015	0.034	0.036	0.039	0.91
ECE-NonP	150	−0.016	0.028	0.030	0.034	0.89
	30	−0.015	0.065	0.067	0.068	0.93
	100	−0.015	0.036	0.036	0.039	0.92
G-estimator	150	−0.015	0.029	0.030	0.034	0.91
	30	−0.015	0.065	0.067	0.068	0.93
	100	−0.015	0.036	0.036	0.039	0.92
EMEE	150	−0.015	0.029	0.030	0.034	0.91
	30	−0.017	0.067	0.067	0.069	0.94
	100	−0.015	0.037	0.036	0.039	0.93
EMEE-NonP	150	−0.016	0.030	0.030	0.034	0.92
	30	−0.003	0.066	0.068	0.068	0.94
	100	−0.002	0.036	0.037	0.037	0.94
DR-EMEE-NonP	150	−0.002	0.030	0.031	0.031	0.95
	30	−0.031	0.062	0.064	0.071	0.92
	100	−0.031	0.034	0.034	0.046	0.85
GEE (ind)	150	−0.032	0.028	0.029	0.043	0.78
	30	−0.031	0.062	0.064	0.071	0.92
	100	−0.031	0.034	0.034	0.046	0.85
GEE (exch)	150	−0.032	0.028	0.029	0.043	0.78

Bias or coverage probability is in bold if it deviates significantly from 0 or 0.95, respectively.

Abbreviations: CP, coverage probability; GEE, generalized estimating equations; RMSE, root mean squared error; SD, standard deviation; SE, standard error.

TABLE 4

Comparison of 8 estimators of the treatment effect moderation |$Z_t$| under Scenario (2).

		|$\beta _0$|					|$\beta _1$|
Estimator	Time length	Bias	SE	SD	RMSE	CP	Bias	SE	SD	RMSE	CP
	30	0.002	0.084	0.088	0.088	0.95	−0.002	0.080	0.081	0.081	0.94
	100	0.000	0.046	0.046	0.046	0.95	0.000	0.044	0.044	0.044	0.95
ECE	150	−0.001	0.038	0.038	0.038	0.95	0.001	0.036	0.035	0.035	0.94
	30	−0.004	0.103	0.089	0.089	0.98	0.000	0.076	0.082	0.082	0.92
	100	−0.005	0.057	0.046	0.047	0.99	0.002	0.041	0.045	0.045	0.94
ECE-NonP	150	−0.005	0.046	0.038	0.038	0.98	0.002	0.034	0.036	0.036	0.93
	30	0.002	0.083	0.086	0.086	0.95	−0.002	0.078	0.080	0.080	0.95
	100	0.000	0.045	0.045	0.045	0.95	0.000	0.043	0.043	0.043	0.95
G-estimator	150	−0.001	0.037	0.037	0.037	0.95	0.001	0.035	0.034	0.034	0.94
	30	0.002	0.083	0.086	0.086	0.95	−0.002	0.078	0.080	0.080	0.95
	100	0.000	0.045	0.045	0.045	0.95	0.000	0.043	0.043	0.043	0.95
EMEE	150	−0.001	0.037	0.037	0.037	0.95	0.001	0.035	0.034	0.034	0.94
	30	−0.004	0.083	0.088	0.088	0.94	0.000	0.076	0.081	0.081	0.94
	100	−0.005	0.045	0.046	0.046	0.95	0.002	0.042	0.044	0.044	0.94
EMEE-NonP	150	−0.005	0.037	0.038	0.038	0.94	0.002	0.034	0.035	0.035	0.94
	30	−0.003	0.086	0.088	0.088	0.95	0.000	0.075	0.081	0.081	0.93
	100	−0.004	0.047	0.046	0.046	0.96	0.002	0.041	0.044	0.044	0.94
DR-EMEE-NonP	150	−0.005	0.038	0.038	0.038	0.95	0.002	0.034	0.035	0.035	0.94
	30	0.002	0.085	0.089	0.089	0.94	−0.002	0.080	0.082	0.082	0.94
	100	0.000	0.047	0.046	0.046	0.96	0.000	0.044	0.044	0.044	0.96
GEE (ind)	150	−0.001	0.038	0.038	0.038	0.95	0.001	0.036	0.035	0.035	0.95
	30	0.002	0.085	0.089	0.089	0.94	−0.002	0.080	0.082	0.082	0.94
	100	0.000	0.047	0.046	0.046	0.96	0.000	0.044	0.044	0.044	0.96
GEE (exch)	150	−0.001	0.038	0.038	0.038	0.95	0.001	0.036	0.035	0.035	0.94

Abbreviations: CP, coverage probability; GEE, generalized estimating equations; RMSE, root mean squared error; SD, standard deviation; SE, standard error.

6 APPLICATION TO DRINK LESS DATA

In this section, we analyze the Drink Less trial data to evaluate the efficacy of delivering push notifications on user engagement with the Drink Less app. Due to the form of intervention, participants were available for the intervention at all times, i.e., |$I_t = 1$|⁠. We exclusively apply EMEE, EMEE-NonP, and DR-EMEE-NonP to the Drink Less data, focusing solely on marginal causal excursion effects, for which ECE, ECE-NonP, and G-estimator have been identified as potentially biased. For EMEE, we adopt a working model |$g(H_t)^{\top }\alpha = \alpha _0 + \alpha _1 Z_t$| for the logarithm of the expected outcome under no treatment. In EMEE-NonP and DR-EMEE-NonP, we leverage 2-part generalized additive models to estimate the conditional mean of proximal outcomes. Baseline and time-varying covariates, including age, gender, employment type, AUDIT score, days since download, and the number of screen views yesterday are employed as control variables in |$g(H_t)$|⁠.

6.1 Primary analysis

In the primary analysis, we estimate the marginal excursion effect of push notifications on the number of screen views. We adopt an analysis model with |$S_t = 1$|⁠, given by

$$\begin{eqnarray} \log \frac{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t, 1} \mid H_t, A_t=1\right)\right\rbrace }{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t, 1} \mid H_t, A_t=0\right)\right\rbrace }=\beta _0. \end{eqnarray}$$

Results in Table 5 show all estimators indicate the effect is statistically significant from zero, with EMEE-NonP and DR-EMEE-NonP yielding marginally reduced SE. Subsequently, we compare the marginal excursion effect between the standard notification and the notification from a new message bank using the following analysis model:

$$\begin{eqnarray} \log \frac{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t, 1} \mid H_t, \boldsymbol {A}_t\right)\right\rbrace }{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t, 1} \mid H_t,\boldsymbol {A}_t=0\right)\right\rbrace }=\beta _{0,1}A_{t,1} + \beta _{0,2}A_{t,2}, \end{eqnarray}$$

where |$A_{t,1} = 1$| denotes the standard notification and |$A_{t,2}=1$| the new notification. Findings in Table 5 affirm the efficacy of both notifications, with the standard notification having a higher treatment effect.

TABLE 5

Marginal excursion effects and effect moderation in the Drink Less micro-randomized trial.

	|$\beta _0$|				|$\beta _1$|
Estimator	Estimate	SE	95% CI	P-value	Estimate	SE	95% CI	P-value
Marginal excursion effects of providing notifications
EMEE	1.110	0.123	(0.869, 1.352)	|$\lt $|⁠.001
EMEE-NonP	1.120	0.111	(0.903, 1.337)	|$\lt $|⁠.001
DR-EMEE-NonP	1.120	0.089	(0.946, 1.294)	|$\lt $|⁠.001
Marginal excursion effects of providing standard notifications
EMEE	1.242	0.132	(0.983, 1.500)	|$\lt $|⁠.001
EMEE-NonP	1.238	0.130	(0.982, 1.493)	|$\lt $|⁠.001
DR-EMEE-NonP	1.240	0.100	(1.042, 1.436)	|$\lt $|⁠.001
Marginal excursion effects of providing new notifications
EMEE	0.965	0.138	(0.694, 1.236)	|$\lt $|⁠.001
EMEE-NonP	0.988	0.129	(0.736, 1.240)	|$\lt $|⁠.001
DR-EMEE-NonP	0.991	0.099	(0.798, 1.184)	|$\lt $|⁠.001
Effect moderation of days since download
EMEE	1.352	0.181	(0.997, 1.707)	|$\lt $|⁠.001	−0.019	0.011	(−0.040, 0.003)	.091
EMEE-NonP	1.555	0.185	(1.187, 1.914)	|$\lt $|⁠.001	−0.034	0.013	(−0.059, −0.009)	.001
DR-EMEE-NonP	1.477	0.153	(1.177, 1.777)	|$\lt $|⁠.001	−0.027	0.010	(−0.046, −0.008)	.001
Effect moderation of the number of screen views yesterday
EMEE	1.195	0.164	(0.873, 1.517)	|$\lt $|⁠.001	−0.023	0.024	(−0.071, 0.024)	.329
EMEE-NonP	1.116	0.122	(0.878, 1.355)	|$\lt $|⁠.001	0.000	0.018	(−0.034, 0.035)	.982
DR-EMEE-NonP	1.118	0.096	(0.930, 1.307)	|$\lt $|⁠.001	0.001	0.017	(−0.032, 0.033)	.974

Abbreviations: SE, standard error; 95% CI, 95% confidence interval.

6.2 Secondary analysis

For the secondary analysis, we examine the effect moderation of providing push notifications on user engagement by setting |$S_t$| to variables such as “days since download” and “the number of screen views yesterday.” The analysis model is

$$\begin{eqnarray} \log \frac{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t, 1} \mid H_t, A_t=1\right) \mid S_t\right\rbrace }{\mathbb {E}\left\lbrace \mathbb {E}\left(Y_{t, 1} \mid H_t, A_t=0\right) \mid S_t\right\rbrace }=\beta _0+\beta _1 S_t. \end{eqnarray}$$

Findings in Table 5 show EMEE-NonP and DR-EMEE-NonP identify “days since download” as a significant moderator, overlooked by EMEE. Furthermore, all estimators reveal “number of screen views yesterday” does not significantly influence the treatment effect.

Based on the findings, push notifications significantly affect user engagement, with the standard notification exhibiting a higher effect. Additionally, over prolonged use, users seem to get habituated, resulting in diminishing treatment effects. This highlights the potential utility of periodically refreshing intervention strategies to maintain user engagement.

7 DISCUSSION

This paper revisits causal excursion effects, contrasting 2 excursions from the current treatment protocol into the future, focusing on zero-inflated count proximal outcomes. We extend ECE and EMEE methods from Qian et al. (2021a)—originally developed for binary outcomes—to zero-inflated count outcomes and introduce novel nonparametric estimators for nuisance function estimation. Notably, zero-inflated count outcomes are relatively unexplored in mHealth, with few discussions on extensions to count outcomes lacking theoretical and empirical investigation. We demonstrate the rate double robustness of ECE-NonP and DR-EMEE-NonP estimators, useful for analyzing observational mHealth data. Furthermore, we establish consistency and asymptotic normality for ECE-NonP, EMEE-NonP, and DR-EMEE-NonP under bidirectional asymptotics (Yu et al., 2023), requiring either sample size or number of decision points to go to infinity.

We summarize a few directions for future research. Firstly, incorporating random effects into causal excursion effects could enable person-specific effects for informing decision-making. Secondly, identifying potential moderators during initial research phases could guide the selection of |$S_t$| from |$H_t$|⁠, currently governed by researcher’s judgment. Lastly, while our approach employs nonparametric methods for estimating nuisance functions, with the causal excursion effect model remaining parametric for interpretability of low-dimensional models, future research could extend into nonparametric methods for the causal excursion effect, thereby introducing enhanced flexibility.

ACKNOWLEDGMENTS

We thank Dr. Claire Garnett, Dr. Olga Perski, Dr. Henry W.W. Potts, and Dr. Elizabeth Williamson for their important contributions to the Drink Less MRT. The authors would also like to thank the anonymous referees, the Associate Editor, and the Editor for their conscientious efforts and constructive comments, which improved the quality of this paper.

FUNDING

Xueqing Liu is supported by PhD student scholarship from the Duke-NUS Medical School, Singapore. Lauren Bell is supported by a PhD studentship funded by the MRC Network of Hubs for Trials Methodology Research (MR/L004933/2-R18). Bibhas Chakraborty would like to acknowledge support from the grant MOE-T2EP20122-0013 from the Ministry of Education, Singapore.

CONFLICT OF INTEREST

None declared.

DATA AVAILABILITY

The Drink Less micro-randomized trial data that support the findings of this paper are publicly available at https://osf.io/mtcfa.