# Dynamics of droplets under electrowetting effect with voltages exceeding the
contact angle saturation threshold

Quoc Vo1<EMAIL_ADDRESS>Tuan Tran1<EMAIL_ADDRESS>1School of Mechanical &
Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue,
639798, Singapore

###### Abstract

Electrowetting-on-dielectric (EWOD) is a powerful tool in many droplet-
manipulation applications with a notorious weakness caused by contact-angle
saturation (CAS), a phenomenon limiting the equilibrium contact angle of an
EWOD-actuated droplet at high applied voltage. In this paper, we study the
spreading behaviours of droplets on EWOD substrates with the range of applied
voltage exceeding the saturation limit. We experimentally find that at the
initial stage of spreading, the driving force at the contact line still
follows the Young-Lippmann law even if the applied voltage is higher than the
CAS voltage. We then theoretically establish the relation between the initial
contact-line velocity and the applied voltage using the force balance at the
contact line. We also find that the amplitude of capillary waves on the
droplet surface generated by the contact line’s initial motion increases with
the applied voltage. We provide a working framework utilising EWOD with
voltages beyond CAS by characterising the capillary waves formed on the
droplet surface and their self-similar behaviours. We finally propose a
theoretical model of the wave profiles taking into account the viscous effects
and verify this model experimentally. Our results provide avenues to utilise
the EWOD effect with voltages beyond CAS threshold and have strong bearing on
emerging applications such as digital microfluidic and ink-jet printing.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see
http://journals.cambridge.org/data/relatedlink/jfm-keywords.pdf for the full
list)

##  1 Introduction

Contact angle of an electrically conductive droplet on an electrode covered by
a dielectric layer reduces when a voltage $U$ is applied between the droplet
and the electrode (see Fig. 1a). This phenomenon, known as electrowetting-on-
dielectric (EWOD), is one of the most effective techniques for droplet
manipulation in three dimensional digital microfluidics (Jun Lee et al., 2012;
Hong & Lee, 2015), anti-icing (Mishchenko et al., 2010), self-cleaning
(Blossey, 2003) and control of droplet deposition (Baret & Brinkmann, 2006).
At low applied voltages, the equilibrium contact angle $\theta_{\rm e}$ after
$U$ is applied directly relates to the voltage $U$ and droplet’s initial
contact angle $\theta_{\rm 0}$ via the so-called Young-Lippmann equation
(Mugele & Baret, 2005):

$\cos\theta_{\rm e}-\cos\theta_{\rm
0}=\frac{\epsilon\epsilon_{0}U^{2}}{2\sigma d},$ (1)

where $\epsilon_{0}$ is the permittivity of free space; $\epsilon$ and $d$
are, respectively, the dielectric constant and thickness of the dielectric
coating, $\sigma$ is surface tension of the droplet. When the applied voltage
is sufficiently high, $\theta_{\rm e}$ fails to follow Eq. 1 and is limited by
a phenomenon known as contact angle saturation (CAS) (see Fig. 1b).

While transient behaviours of droplets under electrowetting actuation within
the saturation limit are well-established (Mugele & Baret, 2005; Vo & Tran,
2018; Vo et al., 2018), the dynamics of droplets with applied voltage
exceeding the CAS threshold remains elusive. Consequently, most applications
utilising electrowetting uses applied voltage within the saturation value,
$U_{\rm s}$, a requirement significantly limiting EWOD’s capabilities,
especially in manipulating small and viscous droplets (Vo & Tran, 2019; Fair,
2007). Therefore, it is of practical importance to explore how droplet
spreading dynamics is affected when the applied voltage is higher than the CAS
threshold.

In this paper, we investigate the transient behaviours of droplets when the
applied voltage $U$ is higher than the CAS voltage $U_{\rm s}$. We focus on
small droplets with radius below capillary length and study their spreading
behaviours in silicone oil medium under electrowetting actuation. We reveal a
direct influence of increasing the applied voltage on enhancing the initial
spreading velocity of droplets even when $U\geq U_{\rm s}$. We also observe
strong capillary waves on the droplet-oil interface and analyse their
occurrence with the enhancement of the initial spreading velocity. We then
propose and verify a theoretical model to explain the relation between the
capillary waves under electrowetting actuation and the applied voltage.

Figure 1:  (a) Schematic showing the EWOD setup. (b) Plot showing the
dependence of the equilibrium contact angle $\theta_{\rm e}$ on the applied
voltage $U$. The solid curve represents Eq. 1 without any fitting parameter.
Contact angle saturation is observed for $U\geq U_{\rm s}\approx 110\,$V.

##  2 Experimental setup, materials and method

In our experiment, we prepare test substrates using indium-tin-oxide (ITO)
glass slides, each covered by a layer of fluoropolymer (Teflon, Dupont) having
thickness $d=2.5{\rm\upmu m}$. We dip an $18\,\mu$m-wire electrode made of
tungsten into the droplet and connect it to the positive terminal of a DC
power supply via a solid state relay (SSR) (see Fig. 1a), while the negative
terminal to the ITO layer. The voltage applied to the two electrodes takes the
form of a step function having amplitude $U$ and duration $T_{\rm p}$; the
duration $T_{\rm p}$ is set sufficiently long to ensure that droplets under
the electrowetting effect reach new equilibrium in every actuation. The
amplitude $U$ is varied between $0\,{\rm V}$ to $190\,{\rm V}$. The contact
angle saturation voltage (CAS) is experimentally determined at $U_{\rm
s}=110\,$V by observing the saturation behaviour of the contact angle when $U$
is gradually increased (Fig. 1b). We note that the observed change of the
equilibrium contact angle $\theta_{\rm e}$ follows the Young-Lippmann equation
(Eq. 1) when $U<U_{\rm s}$ (Fig. 1b, solid curve).

We use a 0.125 M sodium chloride aqueous solution to generate droplets. For
each experiment, we immerse a droplet and the substrate in silicone oil having
viscosity $\mu_{\rm o}=1.8\,{\rm mPas}$. The interfacial tension of the
solution in the silicone oil is measured experimentally at $\sigma=37.2\pm
0.5\,{\rm mPa\,s}$. The temperature of the oil pool is kept at $20\pm
0.5\,^{\circ}{\rm C}$ to maintain consistent experimental conditions. The
radius of the droplet is $R=0.5\,{\rm mm}$, well below the capillary length
$\l_{\rm c}=[\sigma/(\rho-\rho_{\rm o})g]^{1/2}=5.5\,$mm. Here, $g=9.781\,{\rm
m\,s^{-1}}$ is the gravitational acceleration, $\rho=1000\,{\rm kg\,m^{-3}}$
and $\rho_{\rm o}=873\,{\rm kg\,m^{-3}}$ the density of the working liquid and
the oil, respectively. We capture the behaviours of droplets using a high
speed camera (Photron, SAX2) typically running at 5000 frame-per-second.

##  3 Results and Discussions

Figure 2:  (a) Snapshots showing the contact-line behaviours of droplets under
different applied voltages $U$ varying from $60\,{\rm V}$ to $150\,$V. The
scale bar represents 0.5 mm. Capillary waves on droplet’s surface are observed
when $U\geq 90\,$V. (b) Plot showing the contact-line velocity $v_{\rm ct}$
vs. time $t$: open markers for those of voltage below the saturation value,
i.e., $U<110\,$V; solid markers for those of voltage above the saturation
value, i.e., $U>110\,$V. (c) Log-log plot showing the initial contact-line
velocity $v_{\rm i}$ vs. the applied voltage $U$. Inset: Linear plot of
$v_{\rm i}$ vs. $U$ for $U<90\,$V, the voltage’s range in which no capillary
wave is observed on the droplet’s surface.

###  3.1 Spreading dynamics by electrowetting beyond contact angle saturation

In Fig. 2a, we show several series of snapshots of droplets shortly after a
voltage $U$ is applied. Generally, we observe that higher applied voltage
causes more violent spreading behaviours due to larger contact-line tension
imbalance (Mugele & Baret, 2005). For instance, at $U=60\,$V (Fig. 2a, first
panel), the droplet gently spreads and maintains its spherical shape at all
time. At $U=90\,$V (Fig. 2a, second panel), faster contact-line motion
generates capillary waves on the liquid-oil interface and consequently causes
substantial deformation to the droplet’s shape. When the voltage increases
above $120\,$V, the capillary waves become so strong that the deformation at
the liquid-oil interface causes small satellite droplets to eject from the
primary droplets (Fig. 2a, two last panels). In such droplet-ejection
instances, the size of the satellite droplets increases when $U$ varies from
$120\,$V to $150\,$V.

Surprisingly, we also observe from Fig. 2a that droplet deformation increases
with the applied voltage even at $U\geq U_{\rm s}=110\,$V. As the deformation
is caused by the capillary waves originated from the contact-line initial
motion, we hypothesise that CAS causes little effect to the early dynamics of
the contact line. In Fig. 2b, we show how the contact-line velocity $v_{\rm
ct}$ depends on time for several values of $U$ from $30$ V to $170$ V. We
observe that when $U<U_{\rm s}$, increasing $U$ shifts the velocity curves
upwards, indicating that both early-time and late-time spreading dynamics are
affected by the applied voltage (Fig. 2b, open markers). However, when
$U>U_{\rm s}$, escalating $U$ only increases the initial values of $v_{\rm
ct}$. Subsequently, $v_{\rm ct}$ of all the values of $U>U_{\rm s}$ converse
to the same curve (Fig. 2b, solid markers). This observation confirms that the
contact-line initial velocity and the resulting capillary waves on the
droplet’s surface are not limited by CAS; CAS only affects the late-time
spreading dynamics of the droplets.

In Fig. 2c, we show the initial contact-line velocity $v_{\rm i}$, i.e.,
$v_{\rm ct}$ measured at $t=0$, vs. the applied voltage $U$. We observe that
although $v_{\rm i}$ increases with $U$ in the whole range of the applied
voltage, the increasing rate of $v_{\rm i}$ is higher for $U<90\,$V compared
to that for $U\geq 90\,$V. Beyond the critical voltage $U_{\rm c}=90\,$V,
determined where the switch in the increasing rate of $v_{\rm i}$ occurs, it
is possible to observe capillary waves on the droplet-oil interface (see Fig.
2a). The occurrence of the capillary waves when the applied voltage is
increased beyond $U_{\rm c}$ suggests that the initial driving force generated
by the electrowetting effect is opposed by the resistive force at the contact
line for $U<U_{\rm c}$, and by the droplet’s inertia for $U\geq U_{\rm c}$. As
a result, the relation between $v_{\rm i}$ and $U$ can be obtained by
balancing the driving force $f_{\rm e}=\epsilon\epsilon_{0}(2d)^{-1}U^{2}$,
resulted from strong localisation of the electrical field at the contact line
(Mugele & Baret, 2005), and either the resitive force $v_{\rm i}\lambda+f$ for
$U<U_{\rm c}$, or inertia $\rho r_{0}v_{\rm i}^{2}$ for $U\geq U_{\rm c}$.
Here, $\lambda$ is the contact-line friction coefficient described by the
geometric mean of the droplet’s viscosity and the oil’s viscosity (Vo & Tran,
2018), and $f$ is the contact-line elasticity force caused by microscopic
defects on the surface (De Gennes, 1985; Joanny & de Gennes, 1984). If we
assume that $f_{\rm e}=\epsilon\epsilon_{0}(2d)^{-1}U^{2}$ is independent from
CAS in the early-time dynamic, we obtain

$\displaystyle v_{\rm i}$
$\displaystyle=\epsilon\epsilon_{0}(2d\lambda)^{-1}U^{2}-f\lambda^{-1}{}$
$\displaystyle{\rm for}\,\,U\leq U_{\rm c}=90\,{\rm V},$ (2) $\displaystyle
v_{\rm i}$ $\displaystyle=(\epsilon\epsilon_{0})^{1/2}(2d\rho
r_{0})^{-1/2}(U-U_{\rm c})+v_{\rm c}{}$ $\displaystyle{\rm for}\,\,U\geq
U_{\rm c}=90\,{\rm V},$ (3)

where $v_{\rm c}=\epsilon\epsilon_{0}(2d\lambda)^{-1}U_{\rm
c}^{2}-f\lambda^{-1}$ is the critical contact-line velocity above which the
capillary wave is observed on the droplet-oil interface. In our experiment,
$v_{\rm c}=0.3\,{\rm ms^{-1}}$. In Fig. 2c, we plot both Eq. 2 (solid lines in
the main figure and the inset) and Eq. 3 (dashed line). We observe an
excellent agreement between the experimental data and both equations in their
respective validity ranges. Best fit of Eq. 2 to the data in Fig. 2c gives
$\lambda=0.075\pm 0.005\,{\rm Pa\,s}$, in good agreement with the previously
reported values for contact-line friction coefficient in the same systems
($\lambda\approx 0.1\pm 0.01\,{\rm Pa\,s}$) (Vo et al., 2018; Vo & Tran,
2018). The fitting also reveals the contact-line elasticity force $f=(2.1\pm
0.2)\times 10^{-3}\,{\rm Nm^{-1}}$, which is of the same order of the
previously reported value, i.e., $(8.6\pm 0.9)\times 10^{-3}\,{\rm Nm^{-1}}$
in Vo & Tran (2019).

###  3.2 Capillary waves

Fast contact line motion generated by the electrowetting effect when $U\geq
U_{\rm c}$ causes capillary waves on the droplet’s surface. A remarkable
feature of the generated capillary waves is that their amplitude increase with
the applied voltage $U$ even when $U$ is higher than the CAS voltage (see Fig.
2). To exploit this feature in practical applications, e.g., manipulating
small and viscous droplets (Vo & Tran, 2019; Fair, 2007), inducing droplet
jumping by modulating actuation time (Wang et al., 2017), and controlling
droplet ejection (Vo & Tran, 2021), we examine the induced capillary waves on
the water-oil interface to reveal the dependence of the wave amplitude on the
applied voltage.

####  3.2.1 Capillary-wave generation

To induce capillary waves on the surface of a droplet, the contact-line
velocity has to overcome the inertial-capillary velocity $(\sigma/\rho
r_{0})^{1/2}$, formulated using the stabilising velocity of a deformed
droplet’s surface to minimise its curvature (Taylor, 1959). In normal wetting
phenomena, the maximum velocity of the contact line that can be theoretically
achieved is $\sigma/\lambda$. As a result, the condition enabling capillary
waves becomes $\sigma/\lambda\geq(\sigma/\rho r_{0})^{1/2}$, or
$\xi=\lambda/(\sigma\rho r_{0})^{1/2}\leq 1$. Here, the so-called damping
coefficient $\xi$ measures the viscous effects in both the droplet and the
outer oil against inertia for capillary flows. We note that $\xi$ is defined
similarly to the Ohnesorge number ${\rm Oh}=\mu/(\sigma\rho r_{0})^{1/2}$,
where $\mu$ is the liquid viscosity. In wetting phenomena driven by the
electrowetting effect, $v_{\rm i}$ is determined by Eq. 2 and Eq. 3, resulting
in another condition for capillary waves to occur: $v_{\rm i}\geq(\sigma/\rho
r_{0})^{1/2}$. Using the experimental parameters in our experiment, the
condition for the contact-line velocity to generate capillary waves is $v_{\rm
i}\geq(\sigma/\rho r_{0})^{1/2}=0.27\,{\rm ms^{-1}}$, consistent with our
experimental results shown in Fig. 2c, i.e., capillary waves are observable
when $v_{\rm i}\geq 0.3\,{\rm ms^{-1}}$. We note that if we define a Weber-
like number ${\rm We}=v_{\rm i}^{2}\rho r_{0}/\sigma$, the condition for
capillary waves to occur then takes the form ${\rm We}\geq 1$.

Figure 3:  (a) Plot showing the displacement $a$ vs. the wave’s position
$\theta_{0}-\phi$ of the capillary waves propagating along the droplet-oil
interface at different time $t$ when $U=150\,$V is applied. (b) Plot showing
$a$ vs. $\theta_{0}-\phi$ of the capillary waves at $t/\tau=1$ and at
different applied voltages varying from $100$ V to $190$ V. Here, $\tau=(\rho
r_{0}^{3}/\sigma)^{1/2}$ is the inertia-capillary time. The maximum
displacement (at point M) increases with the applied voltage regardless of
CAS.

####  3.2.2 Capillary-wave characteristics

To quantify the amplitude and phase of the capillary waves on a droplet having
spherical-cap radius $r_{\rm s}$ and contact angle $\theta_{0}$, we use a
polar coordinate system $(r,\phi)$ with origin $(r=0,\phi=0)$ at a distance
$r_{s}$ below the droplet’s apex and the polar axis vertically upward (Fig.
3a, inset). For droplets in our experiment, $\theta_{0}\approx\pi$ and $r_{\rm
s}\approx r_{0}$. In this coordinate system, the droplet-oil profile, which is
also the capillary waveform, is $(r_{\rm d},\theta_{0}-\phi)$, where
$\theta_{0}$ is the wave’s position of the contact line at $t=0$. The
displacement of the capillary wave is therefore $a=r_{\rm d}-r_{s}$, measured
by the deviation of the droplet-oil interface from the droplet’s initial shape
($r_{d}=r_{s}$). In Fig. 3a, we show $a$ vs. $\theta_{0}-\phi$ at several
instances of time $t$ when $U=150\,$V is applied. When $t$ increases, the
waveform moves from left to right indicating that the wave propagates further
from the contact line and towards the apex of the droplet ($\phi=0$).
Furthermore, we observe that while the wave’s propagating velocity reduces
with time, its amplitude increases suggesting that the wave’s kinetic energy,
originated from the fast initial motion of the contact line, is progressively
converted to surface energy during its propagation.

In Fig. 3b, we show a plot of $a$ vs. $\theta_{0}-\phi$ at $t=\tau$ for
different values of $U$ varying from $100$ V to $190$ V. Here, $\tau=(\rho
r_{0}^{3}/\sigma)^{1/2}$ is the inertial-capillary time. We observe that the
displacement at point M significantly increases with $U$, while the wave’s
position of the point M does not depend on $U$. This observation is examined
more quantitatively in Fig. 4a and b where we respectively show $a_{\rm(M)}$
and $\theta_{0}-\phi_{\rm(M)}$ vs. $t/\tau$ for all the values of $U$ ranging
from 90 V to 190 V. Here, $a_{\rm(M)}$ and $\theta_{0}-\phi_{\rm(M)}$ denote
the displacement and the position of the waves at point M, respectively. On
the one hand, we observe that the slopes of $a_{\rm(M)}$ vs. $t/\tau$ increase
with $U$ showing that the dependence of $a_{\rm(M)}$ on $t/\tau$ is more
strongly at higher applied voltages. On the other hand,
$\theta_{0}-\phi_{\rm(M)}$ is independent from $U$ at any values of $t/\tau$
indicating that the propagating velocity only depends on the hydrodynamical
properties of the droplets, not the contact-line velocity.

Figure 4: (a) Plot showing the displacement of point M, $a_{\rm(M)}$, vs.
$t/\tau$ for $U$ varying from $90$ V to $190$ V. Inset: Log-log plot of the
data in the main plot showing the power law dependence
$a_{\rm(M)}\sim(t/\tau)^{2/3}$. The legends are given in (b). (b) Plot showing
the wave’s position of point M, $\theta_{0}-\phi_{\rm(M)}$, vs. $t/\tau$ for
$U$ varying from $90$ V to $190$ V. Inset: Log-log plot of the data in the
main plot revealing the power law dependence
$\phi_{0}-\phi_{\rm(M)}\sim(t/\tau)^{2/3}$. (c) Plot showing the dimensionless
displacement $A$ vs. the dimensionless wave’s position $\Phi$ at different
values of $t/\tau$ varying from $0.05$ to $1.4$. We observe self-similar
behaviours of the capillary waves for $0.2\leq t/\tau\leq 1$ (solid curves).
Whereas the capillary waves are not self-similar for $t/\tau<0.2$ and
$t/\tau>1$ (dashed curves).

The scaling laws for the displacement and the position of capillary waves on
the surface of a spreading droplet were proposed by Ding et al. (2012) by
adapting Keller & Miksis (1983)’s self-similarity theory for surface-tension
driven flows:

$a=Ar_{0}\left(\frac{t}{\tau}\right)^{2/3}\,\,{\rm
and\,\,}\theta_{0}-\phi=\Phi\left(\frac{t}{\tau}\right)^{2/3},$ (4)

where $A$ and $\Phi$ denote the dimensionless displacement and the
dimensionless position of the waves, respectively. In the insets of Fig. 4a
and b, we show the log-log plots of $a_{\rm(M)}$ vs. $t/\tau$ and
$\theta_{0}-\phi_{\rm(M)}$ vs. $t/\tau$, respectively. We observe that
$a_{\rm(M)}\sim(t/\tau)^{2/3}$ and
$\theta_{0}-\phi_{\rm(M)}\sim(t/\tau)^{2/3}$, implying that Eq. 4 is also
applicable in our system in which the capillary waves are generated by
electrowetting actuation. To further test Eq. 4, in Fig. 4c, we show the
dimensionless waveforms, i.e., $A=ar_{0}^{-1}(t/\tau)^{-2/3}$ vs.
$\Phi=(\theta_{0}-\phi)(t/\tau)^{-2/3}$ when $t/\tau$ varies from $0.05$ to
$1.4$. We observe that all the waveforms for $0.2\leq t/\tau\leq 1$ collapse
onto a single one confirming the self-similar behaviour of the capillary waves
in this time range. However, we notice that capillary waves are not self-
similar for $t/\tau<0.2$ and $t/\tau>1$ (see the dashed curves in Fig. 4c). At
$t/\tau<0.2$, i.e., a short time after wetting is initiated, the waves are
local at the contact line region (Cox, 1986), whereas, at $t/\tau>1$ the
capillary waves reach the apex of the droplet and start interacting with the
waves coming from the other side of the droplet. As a result, at both extremes
the capillary waves on the droplet’s surface are more complicated and cannot
be described by a single self-similar scaling law.

####  3.2.3 Profiles of capillary waves

With the self-similar behaviour of the capillary waves, it is possible to
reduce the wave’s profile $a(\phi,t)$ at any time $t$ by a time-invariant
dimensionless profile $A(\Phi)$. The first solution for $A(\Phi)$ was proposed
by Keller & Miksis (1983) who modified the calculation of Jeffreys & Jeffreys
(1999) for inviscid, self-similar capillary-driven flows in a fluid wedge.
Later, Billingham (1999) theoretically incorporated viscosity effects using
the capillary-viscosity velocity. However, both Keller & Miksis (1983)’s and
Billingham (1999)’s solutions fail to fit our experimental data. While Keller
& Miksis (1983)’s lacks of the viscous effects, the model of Billingham (1999)
is only applicable at a very early time after the wave is generated, and in a
very small region close to the wave’s source, e.g., the time and length scales
respectively are $\approx 0.2\,$ns and $\approx 10\,$nm for water.

Figure 5: (a) Examples of fitting the experimental data of $A$ vs. $\Phi$ at
$t/\tau=0.25$ to Eq. 11. Similar results are obtained for other values of
$t/\tau$ varying from $0.2$ to $1$. (b) Plot showing the fitting parameters
$k_{\phi}$ (on the left axis) and $\Phi_{0}$ (on the right axis) vs. the
applied voltage $U$. Both $k_{\phi}=2.4\pm 0.7$ and $\Phi_{0}=1.32\pm 0.06$
appear to be constants in the voltage’s range. (c) Plot representing Eq. 11
for different applied voltage varying from $90$ V to $190$ V, and with
$k_{\phi}=2.4$ and $\Phi_{0}=1.32$.

We seek an alternative expression of $A(\Phi)$ taking into account the
viscosity effect in our system. The self-similar profiles $A(\Phi)$ (see the
profiles of $0.2\leq t/\tau\leq 1$ in Fig. 4c) can be approximated by
spatially decaying waves driven by interfacial tension and resisted by inertia
and viscous forces. For such waves, we expect that $A(\Phi)$ is the solution
of an one-dimensional time-invariant wave equation with a spatial damping
coefficient $\xi=\lambda/(\sigma\rho r_{0})^{1/2}$:

$\frac{d^{2}A}{d\Phi^{2}}+2k\xi\frac{dA}{d\Phi}+k^{2}A=0,$ (5)

where $k=2\pi/\lambda_{\rm\Phi}=2\pi$ is the wave number,
$\lambda_{\rm\Phi}=1$ is the dimensionless wavelength. We note that the
damping coefficient $\xi$ does not depend on the strength of the driving force
at the contact line. Our discussion is also limited to the case that $\xi<1$
for the possibility of capillary wave generation as discussed in Section.
3.2.1.

If we denote $\Phi_{0}\geq 0$ the dimensionless position of the contact line,
the profile $A$ in Eq. 5 is required to satisfy the following boundary
conditions:

$\displaystyle A\bigg{\rvert}_{\Phi=\Phi_{0}}$ $\displaystyle=A_{\rm
ct}=\frac{a_{\rm ct}}{r_{0}(t/\tau)^{2/3}},$ (6)
$\displaystyle\frac{dA}{d\Phi}\bigg{\rvert}_{\Phi=\Phi_{0}}$
$\displaystyle=\frac{dA_{\rm ct}}{d\Phi_{\rm ct}}=\frac{\dot{a}_{\rm
ct}r_{0}^{-1}-(2/3)a_{\rm ct}r_{0}^{-1}t^{-1}}{-\dot{\phi}_{\rm
ct}-(2/3)(\theta_{0}-\phi_{\rm ct})t^{-1}},$ (7)

where $a_{\rm ct}=r_{\rm ct}(\sin\phi_{\rm ct})^{-1}-r_{\rm 0}$ is the contact
line’s displacement, $\phi_{\rm ct}$ its wave’s position, $r_{\rm ct}$ the
contact radius, $\dot{a}_{\rm ct}$ denotes the time derivatives of $a_{\rm
ct}$ and $\phi_{\rm ct}$ the time derivatives of $\dot{\phi}_{\rm ct}$. As we
consider the capillary waves generated short-time after the electrowetting
effect is activated, we can approximate $a_{\rm ct}\approx 0$, $\phi_{\rm
ct}\approx\theta_{0}$, $a_{\rm ct}t^{-1}\approx\dot{a}_{\rm ct}\approx v_{\rm
i}\sin\theta_{0}$ and $-(\theta_{0}-\phi_{\rm ct})t^{-1}\approx\dot{\phi}_{\rm
ct}=-k_{\phi}r_{0}^{-1}(\sigma/\rho r_{0})^{1/2}$, where $k_{\phi}$ is a
positive constant. As a result, we obtain

$\displaystyle A\bigg{\rvert}_{\Phi=\Phi_{0}}$ $\displaystyle=0,$ (8)
$\displaystyle\frac{dA}{d\Phi}\bigg{\rvert}_{\Phi=\Phi_{0}}$
$\displaystyle=\frac{v_{\rm i}\sin\theta_{0}}{k_{\phi}(\sigma/\rho
r_{0})^{1/2}}=\frac{{\rm We}^{1/2}\sin\theta_{0}}{k_{\phi}}.$ (9)

The general solution of Eq. 5 is

$A={\rm e}^{-k\xi(\Phi-\Phi_{0})}\left[\alpha\cos k_{\rm
d}(\Phi-\Phi_{0})+\beta\sin k_{\rm d}(\Phi-\Phi_{0})\right].$ (10)

Here, $k_{\rm d}=2\pi(1-\xi^{2})^{1/2}$ is the attenuated wave number due to
the damping coefficient $\xi$. Applying Eq. 8 and Eq. 9 to Eq. 10, we obtain
$\alpha=0$ and $\beta={\rm We}^{1/2}\sin\theta_{0}(k_{\rm d}k_{\phi})^{-1}$.
Consequently, we obtain an explicit expression of $A(\Phi)$:

$A=\frac{{\rm We}^{1/2}\sin\theta_{0}}{k_{\phi}k_{\rm d}}\,{\rm
e}^{-k\xi(\Phi-\Phi_{0})}\,\sin[k_{\rm d}(\Phi-\Phi_{0})].$ (11)

To verify Eq. 11, we use the least-mean-square method to fit the experimental
data of $A$ vs. $\Phi$ to Eq. 11 with the fitting parameters $\Phi_{0}$ and
$k_{\phi}$. Examples of the fittings at $t/\tau=0.25$ and for different
applied voltages $U$ are shown in Fig. 5a. Similar results obtained for
different values of $t/\tau$ varying from $0.2$ to $1$ indicate a good
agreement between Eq. 11 and the experimental data. In Fig. 5b, we show the
values of $k_{\phi}$ (left axis) and $\Phi_{0}$ (right axis) that both can be
approximated as constants, i.e., $k_{\phi}=2.4\pm 0.7$ and $\Phi_{0}=1.32\pm
0.06$, for $U$ varying from 90 V to 190 V. All solutions for $90\leq U\leq
190$ V are plotted in Fig. 5c.

##  4 Conclusions

We have shown that the initial contact-line velocity of a droplet spreading on
a solid substrate under the electrowetting effect is not affected by the
contact angle saturation (CAS). We theoretically derived and experimentally
verified the relation between the initial contact-line velocity and the
applied voltage. The magnitude of the capillary waves on the droplet-oil
interface also depends on the applied voltage regardless of CAS. Based on
theories describing the self-similar behaviours of interfacial-tension-driven
flows, we then proposed a mathematical expression and provided experimental
verification for the far-field capillary-wave profiles incorporating viscous
effects. The proposed model strengthens our current understanding on how the
capillary waves generates, propagates and decays. The model is also useful in
understanding other capillary-wave-driven phenomena such as droplet jumping
(Vahabi et al., 2018), droplet/bubble ejection by coalescence (Zhang et al.,
2015; Zhang & Thoroddsen, 2008; Zhang et al., 2009), bubble-busting jet at
liquid-air interface (Gordillo & Rodríguez-Rodríguez, 2019), and bouncing of
impacting droplets on a solid substrate (Richard et al., 2002; Renardy et al.,
2003).

## Acknowledgments

This study is supported by Nanyang Technological University, the Republic of
Singapore’s Ministry of Education (MOE, Grant No. MOE2018-T2-2-113), and the
Agency for Science, Technology and Research (A*STAR, Grant No. 1523700102).

## Declaration of Interests

The authors report no conflict of interest.

## References

  * Baret & Brinkmann (2006) Baret, J.-C. & Brinkmann, M. 2006 Wettability Control of Droplet Deposition and Detachment. Phys. Rev. Lett. 96, 146106.
  * Billingham (1999) Billingham, J. 1999 Surface-tension-driven flow in fat fluid wedges and cones. J. Fluid Mech. 397, 45–71.
  * Blossey (2003) Blossey, R. 2003 Self-cleaning surfaces — virtual realities. Nat. Mater. 2, 301–306.
  * Cox (1986) Cox, R. G. 1986 The dynamics of the spreading of liquids on a solid-surface. Part 1. Viscous-flow. J. Fluid Mech. 168, 169–194.
  * De Gennes (1985) De Gennes, P. G. 1985 Wetting: Statics and dynamics. Rev. Mod. Phys. 57, 827–863.
  * Ding et al. (2012) Ding, H., Li, E. Q., Zhang, F. H., Sui, Y., Spelt, P. D. M. & Thoroddsen, S. T. 2012 Propagation of capillary waves and ejection of small droplets in rapid droplet spreading. J. Fluid Mech. 697, 92–114.
  * Fair (2007) Fair, R. B. 2007 Digital microfluidics: is a true lab-on-a-chip possible? Microfluid. Nanofluidics 3, 245–281.
  * Taylor (1959) Taylor, G. I. 1959 The dynamics of thin sheets of fluid. III. Disintegration of fluid sheets. Proc. R. Soc. London. Ser. A. Math. Phys. Sci. 253 (1274), 313–321.
  * Gordillo & Rodríguez-Rodríguez (2019) Gordillo, J. M. & Rodríguez-Rodríguez, J. 2019 Capillary waves control the ejection of bubble bursting jets. J. Fluid Mech. 867, 556–571.
  * Hong & Lee (2015) Hong, J. & Lee, S. J. 2015 Detaching droplets in immiscible fluids from a solid substrate with the help of electrowetting. Lab Chip 15 (3), 900–907.
  * Jeffreys & Jeffreys (1999) Jeffreys, H. & Jeffreys, B. 1999 Methods of Mathematical Physics, 3rd edition. Cambridge University Press.
  * Joanny & de Gennes (1984) Joanny, J. F. & de Gennes, P. G. 1984 A model for contact angle hysteresis. J. Chem. Phys. 81 (1), 552–562.
  * Jun Lee et al. (2012) Jun Lee, S., Lee, S. & Kang, K. H. 2012 Droplet jumping by electrowetting and its application to the three-dimensional digital microfluidics. Appl. Phys. Lett. 100 (8), 081604.
  * Keller & Miksis (1983) Keller, J. B. & Miksis, M. J. 1983 Surface Tension Driven Flows. SIAM J. Appl. Math. 43 (2), 268–277.
  * Mishchenko et al. (2010) Mishchenko, L., Hatton, B., Bahadur, V., Taylor, J. A., Krupenkin, T. & Aizenberg, J. 2010 Design of Ice-free Nanostructured Surfaces Based on Repulsion of Impacting Water Droplets. ACS Nano 4 (12), 7699–7707.
  * Mugele & Baret (2005) Mugele, F. & Baret, J.-C. 2005 Electrowetting: from basics to applications. J. Phys. Condens. Matter 17 (28), R705–R774.
  * Renardy et al. (2003) Renardy, Y., Popinet, S., Duchemin, L., Renardy, M., Zaleski, S., Josserand, C., Drumright-Clarke, M. A., Richard, D., Clanet, C. & Quéré, D. 2003 Pyramidal and toroidal water drops after impact on a solid surface. J. Fluid Mech. 484 (484), 69–83.
  * Richard et al. (2002) Richard, D., Clanet, C. & Quéré, D. 2002 Contact time of a bouncing drop. Nature 417 (6891), 811.
  * Vahabi et al. (2018) Vahabi, H., Wang, W., Mabry, J. M. & Kota, A. K. 2018 Coalescence-induced jumping of droplets on superomniphobic surfaces with macrotexture. Sci. Adv. 4 (11), 1–8.
  * Vo et al. (2018) Vo, Q., Su, H. & Tran, T. 2018 Universal Transient Dynamics of Electrowetting Droplets. Sci. Rep. 8 (1), 836.
  * Vo & Tran (2018) Vo, Q. & Tran, T. 2018 Contact line friction of electrowetting actuated viscous droplets. Phys. Rev. E 97 (6), 063101.
  * Vo & Tran (2019) Vo, Q. & Tran, T. 2019 Critical conditions for jumping droplets. Phys. Rev. Lett. 123 (2), 024502.
  * Vo & Tran (2021) Vo, Q. & Tran, T. 2021 Droplet ejection by electrowetting actuation. Appl. Phys. Lett. 118 (2), 161603.
  * Zhang et al. (2009) Zhang, F. H., Li, E. Q. & Thoroddsen, S. T. 2009 Satellite formation during coalescence of unequal size drops. Phys. Rev. Lett. 102 (10), 104502.
  * Zhang et al. (2015) Zhang, F. H., Thoraval, M. J., Thoroddsen, S. T. & Taborek, P. 2015 Partial coalescence from bubbles to drops. J. Fluid Mech. 782, 209–239.
  * Zhang & Thoroddsen (2008) Zhang, F. H. & Thoroddsen, S. T. 2008 Satellite generation during bubble coalescence. Phys. Fluids 20 (2), 022104.
  * Wang et al. (2017) Wang, Z., Ende, D. van den, Pit, A., Lagraauw, R., Wijnperle, D. & Mugele, F. 2017 Jumping drops on hydrophobic surfaces, controlling energy transfer by timed electric actuation. Soft Matter 13, 4856-4863.


# Super-Exponential Regret for UCT, AlphaGo and Variants

Laurent Orseau, Rémi Munos

(Google DeepMind
<EMAIL_ADDRESS>
###### Abstract

We improve the proofs of the lower bounds of Coquelin and Munos (2007) that
demonstrate that UCT can have $\exp(\dots\exp(1)\dots)$ regret (with
$\Omega(D)$ exp terms) on the $D$-chain environment, and that a ‘polynomial’
UCT variant has $\exp_{2}(\exp_{2}(D-O(\log D)))$ regret on the same
environment — the original proofs contain an oversight for rewards bounded in
$[0,1]$, which we fix in the present draft. We also adapt the proofs to
AlphaGo’s MCTS and its descendants (_e.g._ , AlphaZero, Leela Zero) to also
show $\exp_{2}(\exp_{2}(D-O(\log D)))$ regret.

##  1 Introduction

First we present the $D$-chain environment, and then provide lowers bounds for
Polynomial UCT, AlphaZero, and UCT.

##  2 The $D$-chain environment

The $D$-chain environment (Coquelin and Munos, 2007) is as follows — see Fig.
1. Consider a binary tree of depth $D$ with two actions, 1 and 2. After taking
$d<D$ times the action 1 from the root, with 0 reward, taking the action 2
leads to a terminal state with reward $1-(d+1)/D$. After taking $D$ times the
action 1, the next state is a terminal state with reward 1 — this is the
optimal trajectory.

We call $n_{d}$ the node reached after taking $d$ times the action 1, and we
call $n_{d^{\prime}}$ the node reached after taking $d-1$ times the action 1,
then the action 2. The node $n_{0}$ is the root.

12…121121212(D-1)/D(D-2)/D2/D1/D01 Figure 1: The $D$-chain environment. Edge
labels are actions, and node labels are rewards.

######  Remark 1.

Since the environment is deterministic, algorithms could take advantage of
this by not visiting terminal states more than once, meaning that the lower
bounds would not apply to such algorithms. To make matters more interesting,
we can instead assume that there is no terminal state and that the tree is an
infinite binary tree. This forces search algorithms to visit the same nodes
multiple times in search of better rewards elsewhere in the tree. (Making the
environment stochastic would substantially complicate the analysis.) ∎

##  3 Polynomial UCT lower bound

A trajectory is a sequence of states/actions starting from the root and ending
in a terminal state. Let $m_{i,t}$ be the number of trajectories going through
node $n_{i}$, after observing $t$ trajectories. Note that at the root $n_{0}$
we have $m_{0,t}=t$. Define $X_{i,t}$ to be the empirical mean of the rewards
obtained on the $m_{i,t}$ trajectories going through node $n_{i}$.

Coquelin and Munos (2007) propose a variant of UCT (Kocsis and Szepesvári,
2006) as follows. Define, at trajectory $t+1$, for $i\in\\{d,d^{\prime}\\}$
for each $d<D$,

$\displaystyle B_{i,t+1}=X_{i,t}+\sqrt{\frac{\sqrt{m_{d-1,t}}}{m_{i,t}}}\,.$
(1)

If at trajectory $t+1$ the node $n_{d-1}$ is visited, then the node $n_{j}$
with $j\in\operatorname*{argmax}_{i\in\\{d,d^{\prime}\\}}B_{i,t+1}$ is
visited, with tie breaking in favour of $d^{\prime}$. The corresponding action
is $1$ if $j=d$, or $2$ otherwise.

Let $T+1$ be the first step at which the node $n_{D}$ (with the maximum reward
1) is reached. Then all the nodes $n_{d}$ and none of the nodes
$n_{d^{\prime}}$ are visited on this trajectory, so for each $d\leq D$ we
necessarily have $B_{d,T}\geq B_{d^{\prime},T}$. Moreover, due to tie-
breaking, all nodes $n_{d^{\prime}}$ for $d\leq D$ have been visited at least
once, and thus $X_{d^{\prime},T}=1-d/D$. Also note that $X_{d,T}\leq
1-(d+1)/D$. Therefore, for all $d\leq D$,

$\displaystyle X_{d,T}+\sqrt{\frac{\sqrt{m_{d-1,T}}}{m_{d,T}}}\geq
X_{d^{\prime},T}+\sqrt{\frac{\sqrt{m_{d-1,T}}}{m_{d^{\prime},T}}}\,.$ (2)

For the rest of the proof, let us make the subscripts $T$ implicit for visual
clarity. It follows that

$\displaystyle 1-\frac{d+1}{D}+\sqrt{\frac{\sqrt{m_{d-1}}}{m_{d}}}\geq
1-\frac{d}{D}+\sqrt{\frac{\sqrt{m_{d-1}}}{m_{d^{\prime}}}}$ (3)

and thus

$\displaystyle\sqrt{\frac{\sqrt{m_{d-1}}}{m_{d}}}\geq\sqrt{\frac{\sqrt{m_{d-1}}}{m_{d^{\prime}}}}+\frac{1}{D}\,.$
(4)

From this, by dropping the term with $m_{d^{\prime}}$, we deduce that for all
$d\leq D$,

$\displaystyle m_{d-1}\geq\left(\frac{m_{d}}{D^{2}}\right)^{2}\,.$ (5)

######  Remark 2.

If we replace the rewards $(D-d)/D$ in the environment with just $D-d$, then
the original claims follow trivially. Indeed Eq. 4 becomes
$m_{d-1}\geq(m_{d})^{2}$, and using 111Observe that $m_{D-1}=1$ only, which is
not sufficient for bootstrapping the sequence. $m_{D-2}=2=2^{2^{0}}$ we deduce
that $m_{0}\geq 2^{2^{D-2}}$. However, it is often assumed that rewards should
be bounded by 1. Since the algorithm is not scale-invariant, we need to adapt
the analysis. ∎

By applying Eq. 5 recursively we obtain for all $d<D$, writing
$\tilde{D}=D^{2}$ for clarity,

$\displaystyle m_{0}$
$\displaystyle\geq\frac{m_{1}^{2}}{\tilde{D}^{2}}=\frac{1}{\tilde{D}^{2^{1}}}m_{1}^{2^{1}}\geq\frac{1}{\tilde{D}^{2^{1}}}\left(\frac{m_{2}^{2}}{\tilde{D}^{2}}\right)^{2^{1}}=\frac{1}{\tilde{D}^{2^{1}+2^{2}}}m_{2}^{2^{2}}\geq\frac{1}{\tilde{D}^{2^{1}+2^{2}}}\left(\frac{m_{3}^{2}}{\tilde{D}^{2}}\right)^{2^{2}}=\frac{1}{\tilde{D}^{2^{1}+2^{2}+2^{3}}}m_{3}^{2^{3}}$
$\displaystyle\geq\dots\geq\frac{1}{\tilde{D}^{2^{1}+2^{2}+2^{3}+\dots
2^{d}}}m_{3}^{2^{d}}\geq\frac{1}{\tilde{D}^{2^{d+1}}}m_{d}^{2^{d}}$
$\displaystyle=\left(\frac{m_{d}}{D^{4}}\right)^{2^{d}}\,.$ (6)

We want to find a value $\hat{d}$ for $d$ such that $m_{\hat{d}}/D^{4}\geq 2$.
To do so, by dropping the $1/D$ term in Eq. 4 and simplifying, 222In case of
iterative expansion of the tree as is often done in practice, we would have
$m_{d-1}=m_{d}+m_{d^{\prime}}+1$ and thus still $m_{d-1}\geq 2m_{d}$. we
deduce that for all $d\leq D$ we have $m_{d}\leq m_{d^{\prime}}$, and thus
$m_{d-1}\geq 2m_{d}$. Applying this relation recursively gives us that

$\displaystyle m_{\hat{d}}\geq 2^{1}m_{\hat{d}+1}\geq
2^{2}m_{\hat{d}+2}\geq\dots\geq 2^{D-1-\hat{d}}m_{D-1}=2^{D-\hat{d}-1}$ (7)

using $m_{D-1,T}=1$. Now, as said above, we want $m_{\hat{d}}/D^{4}\geq 2$,
which is satisfied when $2^{D-\hat{d}-1}\geq 2D^{4}$, which is satisfied for
$\hat{d}=\left\lfloor D-2-4\log_{2}D\right\rfloor$. Plugging $d=\hat{d}$ into
Eq. 6 we obtain (with $\exp_{2}(x)=2^{x}$)

$\displaystyle T=m_{0}\geq\exp_{2}(\exp_{2}(\left\lfloor
D-2-4\log_{2}D\right\rfloor))\,.$ (8)

For $D=25$, this gives $T\geq 2^{1024}\geq 10^{100}$, which is intractable,
while a simple breadth-first search in a full complete binary tree of depth 25
would take only $2^{25}$ search steps — which is well tractable.

##  4 AlphaZero lower bound

In this section we merely adapt the steps for Polynomial UCT to a different
definition fo $B_{i,t}$, and we choose in particular a different $\hat{d}$.

The action selection of the ‘MCTS’ algorithm in AlphaGo (Silver et al., 2016)
and its successors, _e.g._ , AlphaGo Zero (Silver et al., 2017b), AlphaZero
(Silver et al., 2017a) and LeelaChess Zero, 333https://slides.com/crem/lc0#/9
has the following form for the $D$-chain environment, for
$i\in\\{d,d^{\prime}\\}$ for all $d<D$,

$\displaystyle
B_{i,t+1}=Q_{i,t}+c_{\text{puct}}P_{i}\frac{\sqrt{m_{d-1}}}{m_{i}+1}$ (9)

where $c_{\text{puct}}$ is a small positive constant such as 2 or 4, $P_{i}$
is the policy weight of action $i$ such that $P_{d}+P_{d^{\prime}}=1$ for all
$d\leq D$, and $Q_{i,t}$ is the “combined mean action value”. In what follows
we assume that (i) $Q_{i,t}=X_{i,t}$ as defined above, (ii)
$P_{d}=P_{d^{\prime}}=1/2$ and $c_{\text{puct}}P_{i}=c$ for some constant
$c>0$. That is,

$\displaystyle B_{i,t+1}=X_{i,t}+c\frac{\sqrt{m_{d-1}}}{m_{i}+1}\,.$ (10)

Similarly to Eq. 4, on trajectory $T+1$ where the maximum reward is reached,
we can derive:

$\displaystyle\frac{\sqrt{m_{d-1}}}{m_{d}+1}\geq\frac{\sqrt{m_{d-1}}}{m_{d^{\prime}}+1}+\frac{1}{D}\,.$
(11)

Still assuming that ties are broken in favour of $n_{d^{\prime}}$, we deduce
that for all $d\leq D,m_{d}\leq m_{d^{\prime}}$ and thus $m_{d-1}\geq 2m_{d}$.
Since $m_{d-1}=m_{d}+m_{d^{\prime}}+1$ (due to one node expansion per
trajectory) we have $m_{D-1}=2$ and thus for all $d<D$,

$\displaystyle m_{d}\geq 2^{D-d}\,.$ (12)

From Eq. 10, we also deduce that

$\displaystyle m_{d-1}\geq\left(\frac{m_{d}}{cD}\right)^{2}\,.$ (13)

Similarly to Eq. 6 where we now take $\tilde{D}=cD$, we obtain for all $d<D$

$\displaystyle m_{0}\geq\left(\frac{m_{d}}{c^{2}D^{2}}\right)^{2^{d}}\,.$ (14)

Now we want to choose a $\hat{d}$ such that $m_{\hat{d}}/c^{2}D^{2}\geq 2$,
which is satisfied for $2^{D-\hat{d}}\geq 2c^{2}D^{2}$, which is satisfied for
$\hat{d}=\left\lfloor D-1-2\log_{2}(cD)\right\rfloor$. This gives

$\displaystyle T=m_{0}~{}\geq~{}\exp_{2}(\exp_{2}(\left\lfloor
D-1-2\log_{2}(cD)\right\rfloor))\,.$ (15)

For example, for $c=2$ ($c_{\text{puct}}=4$) and $D=20$ this gives $T\geq
2^{2048}\geq 10^{200}$.

##  5 UCT lower bound

UCT uses the following formula, for $i\in\\{d,d^{\prime}\\}$ for all $d<D$:

$\displaystyle B_{i,t+1}=X_{i,t}+\sqrt{\frac{2\ln m_{d-1,t}}{m_{i,t}}}\,.$
(16)

Similarly to the previous sections, on the trajectory $T+1$ we obtain, for all
$d<D$,

$\displaystyle B_{d^{\prime},T+1}$ $\displaystyle\leq B_{d^{\prime},T}\,,$
(17) $\displaystyle 1-\frac{d}{D}+\sqrt{\frac{2\ln
m_{d-1,T}}{m_{d^{\prime},T}}}$ $\displaystyle\leq
1-\frac{d+1}{D}+\sqrt{\frac{2\ln m_{d-1,T}}{m_{d,T}}}\,,$ (18)

from which we deduce (omitting $T$ subscripts):

$\displaystyle\sqrt{\frac{2\ln m_{d-1}}{m_{d}}}$
$\displaystyle\geq\frac{1}{D}\,,$ (19) $\displaystyle m_{d-1}$
$\displaystyle\geq\exp\left(\frac{m_{d}}{2D^{2}}\right)\,.$ (20)

Observe that starting the recursion with $m_{d}\leq 2D^{2}$ gives at most
$m_{d-1}\geq\exp(1)$ which is less than $2D^{2}$ (for $D\geq 2$) and the
recurrence actually converges to a number close to 1. Hence, for the
exponential behaviour to kick in, we need to start with $m_{d}$ large enough,
that is, we need to find some $\hat{d}$ such that

$\displaystyle\exp\left(\frac{m_{\hat{d}}}{2D^{2}}\right)\geq 2m_{\hat{d}}\,.$
(21)

If $m_{\hat{d}}\geq 4D^{3}$ then it can be shown that the relation above is
true for all $D\geq 3$.

From Eq. 18 we deduce that Eq. 7 still holds, which means that
$m_{\hat{d}}\geq 4D^{3}$ is satisfied when $2^{D-\hat{d}-1}\geq 4D^{3}$, which
is satisfied for $\hat{d}=\left\lfloor D-3\log_{2}D-3\right\rfloor$.

We conclude that

$\displaystyle
m_{0}\geq\exp(\dots\exp(\exp(4D^{3}/2D^{2})/2D^{2})\dots/2D^{2})\,,$ (22)

where the number of exp is $\left\lfloor
D-3\log_{2}D-3\right\rfloor=\Omega(D)$. For $D=16$ already, $m_{0}\geq
e^{e^{25}}$.

## References

  * Coquelin and Munos [2007] Pierre-Arnaud Coquelin and Rémi Munos.  Bandit algorithms for tree search.  In Proceedings of the Twenty-Third Conference on Uncertainty in Artificial Intelligence, UAI’07, page 67–74, Arlington, Virginia, USA, 2007\. AUAI Press.
  * Kocsis and Szepesvári [2006] Levente Kocsis and Csaba Szepesvári.  Bandit based monte-carlo planning.  In Johannes Fürnkranz, Tobias Scheffer, and Myra Spiliopoulou, editors, Machine Learning: ECML 2006, pages 282–293, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg.
  * Silver et al. [2016] David Silver, Aja Huang, Christopher J. Maddison, Arthur Guez, Laurent Sifre, George van den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, Sander Dieleman, Dominik Grewe, John Nham, Nal Kalchbrenner, Ilya Sutskever, Timothy Lillicrap, Madeleine Leach, Koray Kavukcuoglu, Thore Graepel, and Demis Hassabis.  Mastering the game of go with deep neural networks and tree search.  Nature, 529:484–503, 2016.
  * Silver et al. [2017a] David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, et al.  Mastering chess and shogi by self-play with a general reinforcement learning algorithm.  arXiv preprint arXiv:1712.01815, 2017.
  * Silver et al. [2017b] David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, Yutian Chen, Timothy Lillicrap, Fan Hui, Laurent Sifre, George van den Driessche, Thore Graepel, and Demis Hassabis.  Mastering the game of go without human knowledge.  Nature, 550(7676):354–359, October 2017.


plain

# A Risk-Ratio-Based Marginal Sensitivity Model for Causal Effects in
Observational Studies

MD ABDUL BASIT  MAHBUB AHM LATIF<EMAIL_ADDRESS><EMAIL_ADDRESS>Institute of Statistical Research and Training,
Univeristy of Dhaka, Bangladesh  ABDUS S. WAHED
<EMAIL_ADDRESS>Department of Biostatistics and Computational
Biology,
University of Rochester, USA

###### Abstract

In observational studies, the identification of causal estimands depends on
the no unmeasured confounding (NUC) assumption. As this assumption is not
testable from observed data, sensitivity analysis plays an important role in
observational studies to investigate the impact of unmeasured confounding on
the causal conclusions. In this paper, we proposed a risk-ratio-based
sensitivity analysis framework by introducing a modified marginal sensitivity
model for observational studies with binary treatments. We further extended
the proposed framework to the multivalued treatment setting. We then showed
how the point estimate intervals and the corresponding percentile bootstrap
confidence intervals can be constructed efficiently under the proposed
framework. Simulation results suggested that the proposed framework of
sensitivity analysis performs well in the presence of adequate overlap among
the treatment groups. Lastly, we demonstrated our proposed sensitivity
analysis framework by estimating the causal effect of maternal education on
female fertility in Bangladesh.

###### keywords:

Causal Inference, Inverse Probability Weighting, Percentile Bootstrap,
Sensitivity Analysis.

##  section Introduction

Randomized controlled trials (RCTs) are considered to be the gold standard for
estimating the causal effect of a treatment on an outcome of interest.
However, researchers are often limited to investigating causal relationships
using only observational data when conducting randomized experiments is not
feasible. In the absence of randomized treatment assignments, the estimation
of causal effects in observational studies is usually based on a set of
identifiability assumptions. One of the most important such assumptions is the
strong ignorability or the no unmeasured confounding (NUC) assumption, which
implies that we can observe all relevant confounders in an observational
study. Since the NUC assumption is essentially unverifiable using observed
data, it is necessary to perform sensitivity analyses that assess the
robustness of the causal conclusions obtained from observational studies in
the presence of unmeasured confounding.

Sensitivity analysis is recognized as an essential step in the process of
causal inference from observational studies. It refers to the approaches that
investigate how causal conclusions are impacted in the presence of unmeasured
confounding in observational studies. Bind & Rubin (2021), in their guidelines
for the statistical reporting of observational studies, listed sensitivity
analysis as one of the five major steps that need to be carried out before
reporting causal conclusions elicited from observational studies. Many
sponsors and regulatory bodies also nowadays require researchers to conduct
some form of sensitivity analysis while drawing causal inferences from
observational studies (for example, see standard MD-4 of the Patient-Centered
Outcome Research Institute (PCORI) methodology standards for handling missing
data at https://tinyurl.com/kcbkvfwv).

Cornfield et al. (1959) conducted the first formal sensitivity analysis to
investigate the impact of unmeasured confounding on the causal relationship
between smoking and lung cancer. However, this initial sensitivity analysis
framework was applicable to only binary outcomes, and it did not account for
the sampling variation. Rosenbaum and colleagues overcame these limitations
and greatly expanded the theory and methods for sensitivity analysis in a
series of pioneering works (Rosenbaum, 1987; Gastwirth et al., 1998; Rosenbaum
et al., 2002). Recently, sensitivity analysis has gained a lot of research
interest including Bonvini et al. (2022); Dorn et al. (2021); Carnegie et al.
(2016); VanderWeele & Ding (2017); Zhao et al. (2019); Kallus et al. (2019);
Yadlowsky et al. (2022), among others. Most of these proposed sensitivity
analysis frameworks deal with only binary treatments.

Zhao et al. (2019) recently proposed a sensitivity analysis framework for
smooth estimators of causal effects, such as the inverse probability weighting
(IPW) estimators. Their proposed marginal sensitivity model is a natural
modification of Rosenbaum’s sensitivity model (Rosenbaum, 2002) for matched
observational studies, and it quantifies the magnitude of unmeasured
confounding by the odds ratio between the conditional probability of being
treated given the measured confounders (observed propensity score) and
conditional probability of being treated given both the measured and
unmeasured confounders (true propensity score). It is well known that risk
ratios are more consistent with the general intuition than odds ratios, and
hence, are easier to interpret. Therefore, it is desirable to develop
sensitivity analysis frameworks that interpret the sensitivity analysis
results using risk ratios instead of odds ratios.

In this paper, we first modified the odds-ratio-based sensitivity analysis
framework of Zhao et al. (2019) to a risk-ratio-based framework. In
particular, we proposed a modified marginal sensitivity model that measures
the violation of the NUC assumption using the risk ratio between the true
propensity score and the observed propensity score (see Section
section.subsection for the exact definition). The proposed modified marginal
sensitivity model introduces a new implicit sensitivity parameter that
restricts the true propensity scores within its valid range for a given
sensitivity model. We then extended the proposed modified marginal sensitivity
model to the multivalued treatment setting based on the work of Basit et al.
(2023). We showed that the estimation of causal effects for binary and
multivalued treatments can be performed efficiently under the proposed
sensitivity analysis framework.

The rest of the paper is structured as follows: Section section describes
necessary notations and briefly reviews the existing odds-ratio-based marginal
sensitivity model; Section section introduces the proposed modified marginal
sensitivity model for the binary treatment setting; Section section extends
the proposed sensitivity model to multivalued treatment settings; Section
section reports results from simulation studies; Section section demonstrates
sensitivity analysis under the proposed framework using a real data
application; Section section provides some concluding remarks.

##  section The Sensitivity Analysis Framework

###  section.subsection The Potential Outcome Framework

Consider an observational study with a binary treatment $A$ (1, if treated,
and 0, if control), a vector of measured confounders
$\boldsymbol{X}\in\mathscr{X}\subset\mathbb{R}^{d}$, and the observed binary
outcome Y, where $(A,\boldsymbol{X},Y)\sim F_{0}$. We observe
$(A_{1},\boldsymbol{X}_{1},Y_{1})$,
$(A_{2},\boldsymbol{X}_{2},Y_{2}),\ldots,(A_{n},\boldsymbol{X}_{n},Y_{n})$,
which denote data points observed from $n$ i.i.d. units from the true data
generating distribution $F_{0}$. Moreover, let $Y_{i}(0)$ and $Y_{i}(1)$
denote the potential outcomes corresponding to treatment levels $0$ and $1$,
respectively, for each unit $i\in[n]$. Under the stable unit treatment value
assumption (SUTVA) (Rubin, 1978), the observed outcomes can be defined in
terms of potential outcomes as
$Y_{i}=Y_{i}(A_{i})=A_{i}Y_{i}(1)+(1-A_{i})Y_{i}(0)$. In order to estimate
causal effects, we make the following identifiability assumptions:

[strong ignorability or NUC] There is no unmeasured confounding (NUC), i.e.,
$(Y(0),Y(1))\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}A\,|\,\boldsymbol{X}.$ In other words, the set
of observed covariates, $\boldsymbol{X}$ includes all common causes of $A$ and
$Y$.

[positivity or overlap] Each unit has a non-zero probability of receiving the
treatment. That is,
$0<\mathbb{P}_{0}(A=1\,|\,X=x)<1,\thickspace\forall\thickspace\boldsymbol{x}\in\boldsymbol{X}$.

In observational studies with a binary treatment, one of the most commonly
used causal estimands of interest is the average treatment effect (ATE),
$\Delta:=\mathbb{E}_{0}[Y(1)]-\mathbb{E}_{0}[Y(0)]$, where $\mathbb{E}_{0}$
indicates that the expectation is taken over the true data generating
distribution $F_{0}$.

If we define the observed propensity score as
$e_{0}(\boldsymbol{X})=\mathbb{P}_{0}(A=1\,|\,\boldsymbol{X})$, then under
Assumptions section.subsection and section.subsection, a consistent estimator
of ATE $\Delta$ based on inverse probability weighting (IPW) is defined as

$\displaystyle\hat{\Delta}_{\mathrm{IPW}}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{A_{i}Y_{i}}{\hat{e}\left(\boldsymbol{X}_{i}\right)}-\frac{\left(1-A_{i}\right)Y_{i}}{1-\hat{e}\left(\boldsymbol{X}_{i}\right)}\right],$

where $\hat{e}(\boldsymbol{X})$ is a sample estimate of
$e_{0}(\boldsymbol{X})$. It is well known that the IPW estimates become
unstable when the estimated propensity scores are close to $0$ or $1$ (Kang et
al., 2007). Therefore, the stabilized IPW (SIPW) estimator of ATE, obtained by
multiplying the inverse probability weights by the probability of receiving
the actual treatment, is frequently used in practice.

###  section.subsection The Odds-Ratio-Based Marginal Sensitivity Model

Zhao et al. (2019) developed a sensitivity analysis framework that can be used
for smooth estimators of causal effects, e.g., the inverse probability
weighting (IPW) estimators of ATE in observational studies with a binary
treatment. Let us consider an unmeasured confounder $U$ that sums up all
unmeasured confounding present in an observational study. Robins (2002)
suggested that the variable $U$ can be considered as any of the potential
outcomes while defining the conditional probabilities of receiving treatments,
i.e., propensity scores. So, let us denote the unobserved or true propensity
score (conditional on both observed and unobserved confounders) by
$e_{a,0}(\boldsymbol{x},y)=\mathbb{P}_{0}\\{A=a\,|\,\boldsymbol{X}=\boldsymbol{x},Y(a)=y\\}$
for $a\in\\{0,1\\}$ and the observed propensity score (conditional on observed
confounders only) by
$e_{a,0}(\boldsymbol{x})=\mathbb{P}_{0}\\{A=a\,|\,\boldsymbol{X}=\boldsymbol{x}\\}$.
Then, the marginal sensitivity model assumes that

$\displaystyle\frac{1}{\Lambda}\leqslant\operatorname{OR}\big{\\{}e_{a}(\boldsymbol{x},y),e_{a,\boldsymbol{\beta}_{0}}(\boldsymbol{x})\big{\\}}\leqslant\Lambda,\thickspace\forall\thickspace\boldsymbol{x}\in\mathscr{X},y\in\mathbb{R},a\in\\{0,1\\},$
(equation)

where $\operatorname{OR}(p_{1},p_{2})=[p_{1}/(1-p_{1})]/[p_{2}/(1-p_{2})]$ is
the odds ratio, $e_{a,\boldsymbol{\beta}_{0}}(\boldsymbol{x})$ is a parametric
for $e_{a,0}(\boldsymbol{x})$ for $a\in\\{0,1\\}$ using observed confounders,
and $\Lambda\geqslant 1$ is a fixed sensitivity parameter that quantifies the
magnitude of unmeasured confounding by the odds ratio of the true propensity
score $e_{a}(\boldsymbol{x},y)$ and the observed propensity score
$e_{\beta_{a,0}}(\boldsymbol{x})$ for $a\in\\{0,1\\}$. The marginal
sensitivity model can also be defined nonparametrically by replacing
$e_{a,\boldsymbol{\beta}_{0}}(\boldsymbol{x})$ with a nonparametric model
$e_{a,0}(\boldsymbol{x})$ in the sensitivity model defined in Expression
(equation).

Under the marginal sensitivity models, Zhao et al. (2019) used a generalized
mini-max inequality and percentile bootstrap to convert the estimation problem
of causal effects under any given sensitivity model to a linear programming
problem (LPP), which can be solved very efficiently.

##  section The proposed risk-ratio-based Framework of Sensitivity Analysis

###  section.subsection The Modified Marginal Sensitivity Model

For observational studies with binary treatments, the marginal sensitivity
model (equation) presented in Section section measures the violation of the
NUC assumption in terms of the odds ratio between the observed and unobserved
propensity score $e_{0}(\boldsymbol{x})$ and $e_{0}(\boldsymbol{x},y)$. Using
the assumptions and notations defined in Section section, a natural
modification of the marginal sensitivity model (equation) that quantifies the
magnitude of the violation of the NUC assumption by the risk ratio between
$e_{a,0}(\boldsymbol{x})$ and $e_{a,0}(\boldsymbol{x},y)$ could be defined as
$e_{a,0}(\boldsymbol{x},y)\in\mathcal{E}_{\beta_{0}}(\Gamma)$, where

$\displaystyle\mathcal{E}_{\beta_{0}}(\Gamma)=\Big{\\{}e_{a}(\boldsymbol{x},y):\frac{1}{\Gamma}\leq\operatorname{RR}\big{\\{}e_{a,\beta_{0}}(\boldsymbol{x}),e_{a}(\boldsymbol{x},y)\big{\\}}\leq\Gamma,\forall\boldsymbol{x}\in\mathscr{X},y\in\mathbb{R},a\in\\{0,1\\}\Big{\\}},$
(equation)

where $\Gamma\geqslant 1$ is the sensitivity parameter and
$\mathrm{RR}\left(p_{1},p_{2}\right)=p_{1}/p_{2}$ is the risk ratio.

Next, let us define

$\displaystyle
k_{a,\beta_{0}}(\boldsymbol{x})=\log\big{(}e_{a,\beta_{0}}(\boldsymbol{x})\big{)},\;\;\;k_{a,0}(\boldsymbol{x},y)=\log\big{(}e_{a,0}(\boldsymbol{x},y)\big{)},\;\;\;\text{
and }\;\;\;k_{a}(\boldsymbol{x},y)=\log\big{(}e_{a}(\boldsymbol{x},y)\big{)},$

and let
$l_{a,\beta_{0}}(\boldsymbol{x},y)=k_{a,\beta_{0}}(\boldsymbol{x})-k_{a,0}(\boldsymbol{x},y)$
be the log-scale difference of the observed and the unobserved propensity
scores. Similarly, for a postulated sensitivity model
$e_{a}(\boldsymbol{x},y)$, we define
$l_{a}(\boldsymbol{x},y)=k_{a,\beta_{0}}(\boldsymbol{x})-k_{a}(\boldsymbol{x},y)$.
However, under the sensitivity model (equation), the unobserved propensity
score $e_{a}(\boldsymbol{x},y)$ does not lie within $(0,1)$. Because, note
that when $l_{a}(\boldsymbol{x},y)\rightarrow\infty$, then
$\log\big{(}e_{a}(\boldsymbol{x},y)\big{)}\rightarrow-\infty$, and hence,
$e_{a}(\boldsymbol{x},y)\rightarrow 0$. But, as
$l_{a}(\boldsymbol{x},y)\rightarrow-\infty$, then
$\log\big{(}e_{a}(\boldsymbol{x},y)\big{)}\rightarrow\infty$, and hence,
$e_{a}(\boldsymbol{x},y)\rightarrow\infty$. That is, under the sensitivity
model (equation), the unobserved or true propensity score
$e(\boldsymbol{x},y)\in(0,\infty)$.

To circumvent this problem, we introduce an additional sensitivity parameter
in the sensitivity model (equation) that restricts $e_{a}(\boldsymbol{x},y)$
to be in its valid range $(0,1)$ as $l_{a}(x,y)\rightarrow-\infty$.

######  Definition section.theorem (modified marginal sensitivity model).

Fix a parameter $\Gamma_{0}\geqslant 1$. Moreover, define
$\Gamma_{1}=\operatorname{max}\big{\\{}e_{a,\beta_{0}}(\boldsymbol{x}),\Gamma_{0}^{-1}\big{\\}}$.
In observational studies with binary treatments, the modified marginal
sensitivity model assumes that
$e_{a,0}(\boldsymbol{x},y)\in\mathcal{E}_{\beta_{0}}(\Gamma_{0},\Gamma_{1})$,
where

$\displaystyle\mathcal{E}_{\beta_{0}}(\Gamma_{0},\Gamma_{1})=\Big{\\{}e_{a}(\boldsymbol{x},y)$
$\displaystyle:\Gamma_{1}\leq\operatorname{RR}\\{e_{a,\beta_{0}}(\boldsymbol{x}),e_{a}(\boldsymbol{x},y)\\}\leq\Gamma_{0},$
$\displaystyle\;\;\;\;\;\;\;\;\forall\,\boldsymbol{x}\in\mathscr{X},y\in\mathbb{R},a\in\\{0,1\\}\Big{\\}},$
(equation)

and $\mathrm{RR}\left(p_{1},p_{2}\right)=p_{1}/p_{2}$ is the risk ratio.

######  Remark section.theorem.

From the Definition section.theorem, we can see that the sensitivity parameter
$\Gamma_{1}$ implicitly depends on the observed propensity score
$e_{a,\beta_{0}}(\boldsymbol{x})$ and the sensitivity parameter $\Gamma_{0}$.
Therefore, we do not need to explicitly specify the value of $\Gamma_{1}$
while conducting a sensitivity analysis. Consequently, the introduction of the
additional sensitivity parameter $\Gamma_{1}$ in the modified marginal
sensitivity model does not necessarily increase the complexity of conducting
the sensitivity analysis. The main purpose of the parameter $\Gamma_{1}$ is to
restrict the unobserved propensity score $e_{a}(\boldsymbol{x},y)$ within its
valid range $(0,1)$ under a specific sensitivity model.

Now, Under the modified marginal sensitivity model (equation), it is easy to
observe that

$-\gamma_{1}\leq l_{a}(\boldsymbol{x},y)\leq\gamma_{0},$

where $\gamma_{0}=\log(\Gamma_{0})$ and
$\gamma_{1}=\log(\Gamma_{1})=\operatorname{max}\big{\\{}\log(e_{0}(\boldsymbol{x})),-\gamma_{0}\big{\\}}$.
Therefore, it can be shown that the sensitivity model (equation) is similar to
assuming $l_{a}\in\mathcal{L}_{\beta_{0}}\left(\gamma_{0},\gamma_{1}\right)$,
where

$\mathcal{L}_{\beta_{0}}\left(\gamma_{0},\gamma_{1}\right)=\left\\{l_{a}:\mathscr{X}\times\mathbb{R}\rightarrow\mathbb{R}\text{
and }\gamma_{1}\leqslant l_{a}\leqslant\gamma_{0},a\in\\{0,1\\}\right\\}.$

For the remainder of the paper, we consider $l_{a}(\boldsymbol{x},y)$ as
sensitivity models for fixed sensitivity parameters $\gamma_{0}\geqslant 0$.

######  Remark section.theorem.

We can also define the modified marginal sensitivity model (equation)
nonparametrically by replacing $e_{a,\beta_{0}}(\boldsymbol{x})$, and
consequently, $k_{a,\beta_{0}}(\boldsymbol{x})$ with corresponding
nonparametric model $e_{a}(\boldsymbol{x})$ and $k_{a}(\boldsymbol{x})$,
respectively. The statistical methods used in the proposed sensitivity
analysis framework can be applied regardless of the choice of parametric or
non-parametric model for $e_{a}(\boldsymbol{x})$, as long as the model is
smooth enough for a valid bootstrap.

###  section.subsection Estimation of ATE under Modified Marginal Sensitivity
models

In order to estimate the ATE under the modified marginal sensitivity model,
under a specific sensitivity model
$l_{a}\in\mathcal{L}_{\beta_{0}}(\gamma_{0},\gamma_{1})$ for $a\in\\{0,1\\}$,
let us define the shifted propensity scores as

$\displaystyle e^{(l_{a})}_{a}(\boldsymbol{x},y)$
$\displaystyle=\bigg{[}\exp\Big{\\{}l_{a}(\boldsymbol{x},y)-k_{a,\beta_{0}}(\boldsymbol{x})\Big{\\}}\bigg{]}^{-1},$

and the shifted estimand of ATE as

$\displaystyle\Delta^{(l_{0},l_{1})}=$
$\displaystyle\Bigg{\\{}\mathbb{E}_{0}\Bigg{[}\frac{A}{e^{(l_{1})}_{1}(\boldsymbol{X},Y)}\Bigg{]}^{-1}\mathbb{E}_{0}\Bigg{[}\frac{AY}{e^{(l_{1})}_{1}(\boldsymbol{X},Y)}\Bigg{]}\Bigg{\\}}$
$\displaystyle-\Bigg{\\{}\mathbb{E}_{0}\Bigg{[}\frac{1-A}{e^{(l_{0})}_{0}(\boldsymbol{X},Y)}\Bigg{]}^{-1}\mathbb{E}_{0}\Bigg{[}\frac{(1-A)Y}{e^{(l_{0})}_{0}(\boldsymbol{X},Y)}\Bigg{]}\Bigg{\\}}.$
(equation)

Note that, for any $l_{a}\in\mathcal{L}_{\beta_{0}}(\gamma_{0},\gamma_{1})$,
we have

$\displaystyle e^{(l_{a})}_{a}(\boldsymbol{x},y)$
$\displaystyle=\big{[}\exp\big{(}k_{a,\beta_{0}}(\boldsymbol{x})-k_{a}(\boldsymbol{x},y)-k_{a,\beta_{0}}(\boldsymbol{x})\big{)}\big{]}^{-1}$
$\displaystyle=\big{[}\exp\big{(}-\log(e_{a}(\boldsymbol{x},y))\big{)}\big{]}^{-1}$
$\displaystyle=\Bigg{[}\frac{1}{e_{a}(\boldsymbol{x},y)}\Bigg{]}^{-1}=e_{a}(\boldsymbol{x},y).$
(equation)

That is, under any given sensitivity model $l_{a}$, our defined shifted
propensity score is equivalent to the true propensity score
$e_{a}(\boldsymbol{x},y)$. We can estimate these shifted propensity scores
with

$\displaystyle\hat{e}^{(l_{a})}_{a}(\boldsymbol{x},y)$
$\displaystyle=\bigg{[}\exp\Big{\\{}l_{a}(\boldsymbol{x},y)-\hat{k}_{a}(\boldsymbol{x})\Big{\\}}\bigg{]}^{-1},$

where
$\hat{k}_{a}(\boldsymbol{x})=\hat{k}_{a,\hat{\boldsymbol{\beta}}}(\boldsymbol{x})=\log\big{(}\hat{e}_{a,\beta}(\boldsymbol{x})\big{)}$
and $\hat{e}_{a,\beta}(\boldsymbol{x})$ is a parametric estimate of
$e_{a,\beta}(\boldsymbol{x})$ for $a\in\\{0,1\\}$.

Consequently, we can define the stabilized IPW (SIPW) estimate of
$\Delta^{(l_{0},l_{1})}$ as

$\displaystyle\hat{\Delta}^{(l_{0},l_{1})}=$
$\displaystyle\Bigg{\\{}\Bigg{[}\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}}{\hat{e}^{(l_{1})}_{1}\left(\boldsymbol{X}_{i},Y_{i}\right)}\Bigg{]}^{-1}\cdot\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}Y_{i}}{\hat{e}^{(l_{1})}_{1}\left(\boldsymbol{X}_{i},Y_{i}\right)}\Bigg{\\}}$
$\displaystyle-\Bigg{\\{}\Bigg{[}\frac{1}{n}\sum_{i=1}^{n}\frac{1-A_{i}}{\hat{e}^{(l_{0})}_{0}\left(\boldsymbol{X}_{i},Y_{i}\right)}\Bigg{]}^{-1}\cdot\frac{1}{n}\sum_{i=1}^{n}\frac{(1-A_{i})Y_{i}}{\hat{e}^{(l_{0})}_{0}\left(\boldsymbol{X}_{i},Y_{i}\right)}\Bigg{\\}}$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}A_{i}Y_{i}\>\hat{e}^{(l_{1})}_{1}\left(\boldsymbol{X}_{i},Y_{i}\right)}{\sum_{i=1}^{n}A_{i}\>\hat{e}^{(l_{1})}_{1}\left(\boldsymbol{X}_{i},Y_{i}\right)}-\frac{\sum_{i=1}^{n}(1-A_{i})Y_{i}\>\hat{e}^{(l_{0})}_{0}\left(\boldsymbol{X}_{i},Y_{i}\right)}{\sum_{i=1}^{n}(1-A_{i})\>\hat{e}^{(l_{0})}_{0}\left(\boldsymbol{X}_{i},Y_{i}\right)}$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}A_{i}Y_{i}\big{[}\exp\\{l_{1}(\boldsymbol{X}_{i},Y_{i})-\hat{k}_{1}(\boldsymbol{X}_{i})\\}\big{]}}{\sum_{i=1}^{n}A_{i}\big{[}\exp\\{l_{1}(\boldsymbol{X}_{i},Y_{i})-\hat{k}_{1}(\boldsymbol{X}_{i})\\}\big{]}}$
$\displaystyle-\frac{\sum_{i=i}^{n}(1-A_{i})Y_{i}\big{[}\exp\\{l_{0}(\boldsymbol{X}_{i},Y_{i})-\hat{k}_{0}(\boldsymbol{X}_{i})\\}\big{]}}{\sum_{i=1}^{n}(1-A_{i})\big{[}\exp\\{l_{0}(\boldsymbol{X}_{i},Y_{i})-\hat{k}_{0}(\boldsymbol{X}_{i})\\}\big{]}}.$
(equation)

Zhao et al. (2019) showed that estimation problems such as the one given in
Equation (equation) can be transformed to a linear fractional programming
(LFP) problem using Charnes-Cooper transformation (Charnes & Cooper, 1962).
Before defining the LFP for Equation (equation), let us simplify the notations
by assuming that the first $m\leqslant n$ units are in the treatment group
$(A=1)$ and the rest of the observations are in the control group $(A=0)$, and
that the outcomes are sorted in decreasing order among the first $m$ units and
the other $n-m$ units. Finally, we can convert Equation (equation) to the
following LFP

min or max
$\displaystyle\;\;\;\;\frac{\sum_{i=1}^{n}Y_{i}\big{[}z_{i}\exp\\{-\hat{k}_{1}(\boldsymbol{X}_{i})\\}\big{]}}{\sum_{i=1}^{m}\big{[}z_{i}\exp\\{-\hat{k}_{1}(\boldsymbol{X}_{i})\\}\big{]}}-\frac{\sum_{i=m+1}^{n}Y_{i}\big{[}z_{i}\exp\\{-\hat{k}_{0}(\boldsymbol{X}_{i})\\}\big{]}}{\sum_{i=m+1}^{n}\big{[}z_{i}\exp\\{-\hat{k}_{0}(\boldsymbol{X}_{i})\\}\big{]}}$
(equation) subject to $\displaystyle\;\;\;\;\Gamma_{1i}\leqslant
z_{i}\leqslant\Gamma_{0},\;\;\;\;\;\;\text{ for }1\leqslant i\leqslant n,$

where
$z_{i}=\exp\big{\\{}l_{1}\left(\boldsymbol{X}_{i},Y_{i}\right)\big{\\}}$, for
$1\leqslant i\leqslant m,$ and
$z_{i}=\exp\big{\\{}l_{0}\left(\boldsymbol{X}_{i},Y_{i}\right)\big{\\}},$ for
$m+1\leqslant i\leqslant n$. The LFP defined in Equation (equation) can
further be transformed to a linear programming problem (LPP), which can be
solved efficiently. The solution of the LPP (equation) yields a partially
identified point estimate interval of ATE under a specific sensitivity model
$l_{a}\in\mathcal{L}_{\beta_{0}}\left(\gamma_{0},\gamma_{1}\right)$ and
$a\in\\{0,1\\}$. We can also obtain a $100(1-\alpha)\%$ asymptotic confidence
interval for ATE under a postulated sensitivity model $l_{a}$ using a
percentile bootstrap approach (Zhao et al., 2019).

##  section Extension to Observational Studies with Multivalued Treatments

In this section, we extend our proposed sensitivity analysis framework to the
multivalued treatment setting. Suppose we observe i.i.d.
$(A_{i},\boldsymbol{X}_{i},Y_{i})_{i=1}^{n}$ from an observational study with
$J>2$ treatment levels, where $A_{i}\in\mathcal{A}=\\{1,2,\ldots,J\\}$ with
the corresponding set of potential outcomes
$\mathcal{Y}_{i}=\\{Y_{i}(1),Y_{i}(2),\ldots,Y_{i}(J)\\}$, and
$\boldsymbol{X}_{i}\in\mathscr{X}\subset\mathbb{R}^{d}$ is a vector of
observed confounders for each subject $i\in[n]$. Let us also define a
treatment indicator $D_{i}(a)$ ($1$ if $A_{i}=a$, $0$ otherwise) for
$a\in\mathcal{A}$.

The identifiability assumptions for observational studies with multivalued
treatments are almost equivalent to those in the binary treatment setting. We
assume the overlap assumption that implies
$\mathbb{P}_{0}(A=a|\boldsymbol{X}=\boldsymbol{x})>0$ for all
$a\in\mathcal{A}$. However, instead of the strong ignorability assumption, we
assume weak ignorability (Imbens, 2000).

[weak ignorability] $D(a)\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}Y(a)\thinspace|\thinspace\boldsymbol{X},\thickspace\forall\thinspace
a\in\mathcal{A}.$

Imbens (2000) extended the concept of propensity scores from binary treatments
to multivalued treatments introducing the generalized propensity score (GPS),
which is defined as
$r_{a,0}(\boldsymbol{x})=\mathbb{P}_{0}(A=a|\boldsymbol{X}=\boldsymbol{x})$
for $a\in\mathcal{A}$. As in Section section, let us further denote the
unobserved or true propensity score as
$r_{a,0}(\boldsymbol{x},y)=\mathbb{P}_{0}(A=a|\boldsymbol{X}=\boldsymbol{x},Y(a)=y)$.

Based on these assumptions and notations, Basit et al. (2023) recently
proposed a sensitivity analysis framework for the multivalued treatment
setting extending the framework of Zhao et al. (2019) for binary treatments.
Using similar ideas, we propose the following risk-ratio-based modified
marginal sensitive model for observational studies with multivalued
treatments.

######  Definition section.theorem (modified marginal sensitivity model for
multivalued treatments).

For fixed $\Gamma_{0}\geqslant 1$ and
$\Gamma_{1}=\operatorname{max}\big{\\{}r_{a,\beta_{0}}(\boldsymbol{x}),\Gamma_{0}^{-1}\big{\\}}$,
the modified marginal sensitivity model for multivalued treatments assumes
that

$\displaystyle\Gamma_{1}\leqslant\mathrm{RR}\big{\\{}r_{a,\beta_{0}}(\boldsymbol{x}),r_{a}(\boldsymbol{x},y)\big{\\}}\leqslant\Gamma_{0},\thickspace\forall\thickspace\boldsymbol{x}\in\mathscr{X},a\in\mathcal{A},y\in\mathbb{R},$
(equation)

where $\Gamma_{0}$ and $\Gamma_{1}$ are sensitivity parameters and
$\mathrm{RR}\left(p_{1},p_{2}\right)=p_{1}/p_{2}$ is the risk ratio.

Next, as in section section.subsection, let us define

$\displaystyle
k_{a,\beta_{0}}(\boldsymbol{x})=\log\big{(}r_{a,\beta_{0}}(\boldsymbol{x})\big{)}\;\;\;\text{
and }\;\;\;k_{a}(\boldsymbol{x},y)=\log\big{(}r_{a}(\boldsymbol{x},y)\big{)},$

and let
$l_{a}(\boldsymbol{x},y)=k_{a,\beta_{0}}(\boldsymbol{x})-l_{a}(\boldsymbol{x},y)$
for all $a\in\mathcal{A}$ be the log-scale difference between the observed and
the unobserved GPSs under any specified sensitivity model
$r_{a}(\boldsymbol{x},y)$. Then, it can be shown that the sensitivity model
(equation) is equivalent to assuming
$l_{a}\in\mathcal{L}_{\beta_{0}}(\gamma_{0},\gamma_{1})$, where

$\displaystyle\mathcal{L}_{\beta_{0}}(\gamma_{0},\gamma_{1})$
$\displaystyle=\big{\\{}l_{a}:\mathscr{X}\times\mathbb{R}\rightarrow\mathbb{R}\text{
and }\gamma_{1}\leqslant l_{a}\leqslant\gamma_{0},\;a\in\mathcal{A}\big{\\}},$

$\gamma_{0}=\log(\Gamma_{0})$, and
$\gamma_{1}=\log(\Gamma_{1})=\operatorname{max}\big{(}\log(r_{a,\beta_{0}}(\boldsymbol{x})),-\gamma_{0}\big{)}$
are the sensitivity parameters.

In order to define estimands of causal effects in the multivalued treatment
setting, we use a general class of additive causal estimands proposed by Basit
et al. (2023). This class of estimands is based on inverse probability
weighting and is defined as

$\tau(\boldsymbol{c})=\sum_{a=1}^{J}c_{a}m(a)=\sum_{a=1}^{J}c_{a}\mathbb{E}_{0}\Bigg{[}\frac{Y\thinspace
D(a)}{r_{a}(\boldsymbol{X})}\Bigg{]},$ (equation)

where $\boldsymbol{c}=(c_{1},\ldots,c_{J})^{\prime}$ is a vector of contrasts,
$m(a)=\mathbb{E}_{0}[Y(a)]$ is the average potential outcome for
$a\in\mathcal{A}$. These estimands encompass many commonly used causal
estimands of interest for multivalued treatments, such as the pairwise average
treatment effects (PATEs).

Based on the identifiability assumptions defined for multivalued treatments,
we define the shifted estimand of the average potential outcome $m(a)$ under a
specified sensitivity model $l_{a}$ as

$\displaystyle
m^{(h_{a})}(a)=\mathbb{E}_{0}\Bigg{[}\frac{YD(a)}{r_{a}(\boldsymbol{X})}\Bigg{]}=\mathbb{E}_{0}\Bigg{[}\frac{YD(a)}{r_{a}(\boldsymbol{X},Y)}\Bigg{]}=\mathbb{E}_{0}\Bigg{[}\frac{YD(a)}{r^{(h_{a})}(\boldsymbol{X},Y)}\Bigg{]},$
(equation)

where
$r^{(l_{a})}(\boldsymbol{x},y)=\big{[}\exp\big{(}l_{a}(\boldsymbol{x},y)-k_{a,\beta_{0}}(\boldsymbol{x})\big{)}\big{]}^{-1}$
is the shifted GPS. Consequently, under the modified marginal sensitivity
model (equation), the shifted causal estimand becomes

$\displaystyle\tau^{(l_{a})}(\boldsymbol{c})=\sum_{a=1}^{J}c_{a}m^{(l_{a})}(a).$
(equation)

Similar to the estimation of the shifted ATE defined in Equation (equation)
for binary treatments, we can show that the estimation of the causal estimand
can be converted to a linear programming problem (LPP) that allows us to
efficiently estimate the partially identified point estimate intervals of
$\tau(\boldsymbol{c})$ for any
$l_{a}\in\mathcal{L}_{\beta_{0}}(\gamma_{0},\gamma_{1})$. The percentile
bootstrap approach of Zhao et al. (2019) is also applicable for the
computation of $100(1-\alpha)\%$ asymptotic confidence intervals for
$\tau(\boldsymbol{c})$ (Basit et al., 2023).

##  section Simulation Study

We conducted simulated studies to investigate the performance of the proposed
sensitivity analysis framework in the multivalued treatment setting. Our data
generating mechanism and simulation settings are equivalent to Basit et al.
(2023). We simulate three covariates as
$X_{i1}\sim\operatorname{Bernoulli}(0.5)$, $X_{i2}\sim\operatorname{U}(-1,1)$,
and $X_{i3}\sim\operatorname{N}(0,0.5)$. For each $i\in[n]$, the covariate
vector then becomes $\boldsymbol{X}_{i}=(1,X_{i1},X_{i2},X_{i3})^{T}$. We
simulated the treatment assignment mechanism using the following multinomial
distribution

$\big{(}D_{i}(1),D_{i}(2),D_{i}(3)\big{)}\;\big{|}\;\boldsymbol{X}_{i}\sim\operatorname{Multinom}\big{(}r_{1}(\boldsymbol{X}_{i}),r_{2}(\boldsymbol{X}_{i}),r_{3}(\boldsymbol{X}_{i})\big{)},$

where $D_{i}(a)$ is the treatment indicator and

$r_{a}(\boldsymbol{X}_{i},Y_{i})=r_{a}(\boldsymbol{X}_{i})=\frac{\exp\big{(}\boldsymbol{X}_{i}^{T}\beta_{a}\big{)}}{\sum_{a^{\prime}=1}^{3}\exp\big{(}\boldsymbol{X}_{i}^{T}\beta_{a^{\prime}}\big{)}}$

denotes the complete or unobserved GPSs for $a\in\\{1,2,3\\}$ with
$\beta_{1}=(0,0,0,0)^{T}$, $\beta_{2}=k_{2}\times(0,1,1,1)^{T}$, and
$\beta_{3}=k_{3}\times(0,1,1,-1)^{T}$. In order to assess the influence of the
degree of overlap on the point estimates and confidence intervals of our
proposed frameworks, we considered two simulation scenarios. In the first
scenario, we set $(k_{2},k_{3})=(0.1,-0.1)$ to simulate a scenario with
adequate overlap in the covariates, and in the second scenario, we set
$(k_{2},k_{3})=(3,3)$ to induce lack of overlap. The potential outcomes are
generated from the following multinomial distribution

$\big{(}Y_{i}(1),Y_{i}(2),Y_{i}(3)\big{)}\;\big{|}\;\boldsymbol{X}_{i}\sim\operatorname{Multinom}\big{(}p_{Y_{1}}(\boldsymbol{X}_{i}),p_{Y_{2}}(\boldsymbol{X}_{i}),p_{Y_{3}}(\boldsymbol{X}_{i})\big{)},$

where $Y_{i}(a)$ is the potential outcome for treatment level $a$ and

$p_{Y_{a}}(\boldsymbol{X}_{i})=\mathbb{P}(Y(a)=1\big{|}\boldsymbol{X}_{i})=\frac{\exp\big{(}\boldsymbol{X}_{i}^{T}\delta_{a}\big{)}}{\sum_{a^{\prime}=1}^{3}\exp\big{(}\boldsymbol{X}_{i}^{T}\delta_{a^{\prime}}\big{)}}$

with $\delta_{1}=(1,1,1,1)^{T}$, $\delta_{2}=(1,1,-1,1)^{T}$, and
$\delta_{3}=(1,1,1,-1)^{T}$. The observed outcome was then obtained as
$Y_{i}=\sum_{a=1}^{J}D_{i}(a)Y_{i}(a)$ for each subject $i\in[n]$. We
simulated $1000$ datasets with sample size $n=750$ for each scenario and
estimate the interval of point estimates and confidence intervals for the
pairwise ATEs using our proposed sensitivity analysis framework. The pairwise
ATEs are denoted by $\tau_{i,j}$, for ${i,j}\in\\{1,2,3\\}$. We considered six
values of the sensitivity parameters $\gamma_{0}$, namely,
$\gamma_{0}=\\{0,0.1,0.2,0.5,1,2\\}$. The true partially identified intervals
were obtained under each simulation scenario using large scale numerical
approximations. Since the treatment under consideration is multivalued with
three treatment levels, the observed generalized propensity scores (GPSs) are
modeled using the multinomial logit regression model.

Table table: Simulation results for the sensitivity analysis of pairwise ATEs
under the RR framework in an observational study with three treatment levels

| % Bias |
---|---|---
Scenario | ATE | $\gamma_{0}$ | $\Gamma_{0}$ | Lower | Upper |  Non- coverage |  Partially identified interval |  Point estimate interval |  Confidence interval
|  | 0.0 | 1.00 | $2.79$ | $2.79$ | 0.103 | $(-0.050,-0.050)$ | $(-0.048,-0.048)$ | $(-0.115,0.017)$
|  | 0.1 | 1.11 | $1.76$ | $2.88$ | 0.099 | $(-0.137,0.038)$ | $(-0.136,0.040)$ | $(-0.201,0.106)$
|  | 0.2 | 1.22 | $2.64$ | $1.72$ | 0.108 | $(-0.223,0.126)$ | $(-0.223,0.127)$ | $(-0.285,0.193)$
|  | 0.5 | 1.65 | $4.79$ | $0.76$ | 0.116 | $(-0.460,0.375)$ | $(-0.459,0.376)$ | $(-0.512,0.436)$
|  | 1.0 | 2.72 | $11.14$ | $-6.18$ | 0.123 | $(-0.744,0.694)$ | $(-0.741,0.693)$ | $(-0.771,0.729)$
| $\tau_{1,2}$ | 2.0 | 7.39 | $5.03$ | $-7.64$ | 0.095 | $(-0.900,0.882)$ | $(-0.900,0.882)$ | $(-0.915,0.900)$
|  | 0.0 | 1.00 | $1.49$ | $1.49$ | 0.108 | $(-0.034,-0.034)$ | $(-0.035,-0.035)$ | $(-0.102,0.032)$
|  | 0.1 | 1.11 | $2.56$ | $1.42$ | 0.107 | $(-0.121,0.053)$ | $(-0.121,0.053)$ | $(-0.186,0.119)$
|  | 0.2 | 1.22 | $4.83$ | $2.18$ | 0.109 | $(-0.207,0.139)$ | $(-0.207,0.140)$ | $(-0.271,0.205)$
|  | 0.5 | 1.65 | $6.45$ | $-1.08$ | 0.113 | $(-0.444,0.385)$ | $(-0.443,0.385)$ | $(-0.499,0.444)$
|  | 1.0 | 2.72 | $9.69$ | $-5.88$ | 0.117 | $(-0.735,0.699)$ | $(-0.733,0.697)$ | $(-0.766,0.734)$
| $\tau_{1,3}$ | 2.0 | 7.39 | $10.19$ | $-1.84$ | 0.102 | $(-0.903,0.885)$ | $(-0.903,0.885)$ | $(-0.918,0.904)$
|  | 0.0 | 1.00 | $1.22$ | $1.22$ | 0.097 | $(0.015,0.015)$ | $(0.016,0.016)$ | $(-0.052,0.084)$
|  | 0.1 | 1.11 | $0.78$ | $-1.32$ | 0.098 | $(-0.075,0.106)$ | $(-0.074,0.106)$ | $(-0.142,0.173)$
|  | 0.2 | 1.22 | $2.95$ | $-0.14$ | 0.105 | $(-0.165,0.194)$ | $(-0.164,0.195)$ | $(-0.230,0.260)$
|  | 0.5 | 1.65 | $2.88$ | $-2.97$ | 0.111 | $(-0.414,0.440)$ | $(-0.413,0.440)$ | $(-0.471,0.495)$
|  | 1.0 | 2.72 | $8.49$ | $-5.93$ | 0.121 | $(-0.722,0.735)$ | $(-0.719,0.733)$ | $(-0.753,0.763)$
I | $\tau_{2,3}$ | 2.0 | 7.39 | $4.67$ | $-1.94$ | 0.085 | $(-0.897,0.898)$ | $(-0.897,0.899)$ | $(-0.913,0.914)$
|  | 0.0 | 1.00 | $-13.40$ | $-13.40$ | 0.188 | $(-0.047,-0.047)$ | $(-0.084,-0.084)$ | $(-0.261,0.097)$
|  | 0.1 | 1.11 | $-10.14$ | $-16.80$ | 0.197 | $(-0.133,0.040)$ | $(-0.166,0.000)$ | $(-0.332,0.184)$
|  | 0.2 | 1.22 | $-7.29$ | $-20.24$ | 0.204 | $(-0.214,0.124)$ | $(-0.242,0.080)$ | $(-0.397,0.263)$
|  | 0.5 | 1.65 | $0.97$ | $-29.69$ | 0.239 | $(-0.427,0.352)$ | $(-0.442,0.300)$ | $(-0.561,0.475)$
|  | 1.0 | 2.72 | $9.99$ | $-42.90$ | 0.287 | $(-0.673,0.638)$ | $(-0.676,0.591)$ | $(-0.747,0.713)$
| $\tau_{1,2}$ | 2.0 | 7.39 | $16.52$ | $-54.85$ | 0.393 | $(-0.891,0.893)$ | $(-0.888,0.872)$ | $(-0.912,0.911)$
|  | 0.0 | 1.00 | $-14.93$ | $-14.93$ | 0.205 | $(-0.031,-0.031)$ | $(-0.072,-0.072)$ | $(-0.244,0.102)$
|  | 0.1 | 1.11 | $-12.24$ | $-17.45$ | 0.207 | $(-0.117,0.052)$ | $(-0.155,0.007)$ | $(-0.314,0.186)$
|  | 0.2 | 1.22 | $-8.85$ | $-20.76$ | 0.219 | $(-0.201,0.132)$ | $(-0.234,0.083)$ | $(-0.379,0.265)$
|  | 0.5 | 1.65 | $-2.15$ | $-29.86$ | 0.240 | $(-0.422,0.349)$ | $(-0.443,0.295)$ | $(-0.556,0.463)$
|  | 1.0 | 2.72 | $7.20$ | $-42.13$ | 0.269 | $(-0.683,0.626)$ | $(-0.689,0.578)$ | $(-0.755,0.700)$
| $\tau_{1,3}$ | 2.0 | 7.39 | $14.25$ | $-54.48$ | 0.371 | $(-0.903,0.885)$ | $(-0.902,0.863)$ | $(-0.921,0.902)$
|  | 0.0 | 1.00 | $-1.32$ | $-1.32$ | 0.128 | $(0.015,0.015)$ | $(0.014,0.014)$ | $(-0.108,0.134)$
|  | 0.1 | 1.11 | $-1.31$ | $-2.93$ | 0.122 | $(-0.072,0.100)$ | $(-0.074,0.098)$ | $(-0.195,0.218)$
|  | 0.2 | 1.22 | $0.49$ | $-3.65$ | 0.131 | $(-0.156,0.177)$ | $(-0.158,0.176)$ | $(-0.274,0.290)$
|  | 0.5 | 1.65 | $4.48$ | $-8.21$ | 0.136 | $(-0.378,0.379)$ | $(-0.380,0.375)$ | $(-0.476,0.472)$
|  | 1.0 | 2.72 | $10.47$ | $-14.50$ | 0.171 | $(-0.643,0.622)$ | $(-0.642,0.616)$ | $(-0.702,0.679)$
II | $\tau_{2,3}$ | 2.0 | 7.39 | $18.63$ | $-20.99$ | 0.225 | $(-0.882,0.862)$ | $(-0.879,0.857)$ | $(-0.899,0.881)$

  * †

The first four columns represent the simulation settings: Scenario and ATE are
the simulation scenarios and pairwise ATEs of interest, respectively;
$\lambda$ is the sensitivity parameter and $\Lambda=exp(\lambda)$. The next
five columns are the results: respectively percentage average bias (in SD
units) in the lower and upper bound of the point estimate interval; the non-
coverage rate of the confidence intervals under the proposed framework (the
desired non-coverage rate is $10\%$); the partially identified point estimate
intervals of the pairwise ATEs; the median interval for the point estimates;
the median confidence interval.

We simulate $1000$ datasets for each scenario and estimate the interval of
point estimates and construct $90\%$ confidence intervals for the pairwise
ATEs under the proposed sensitivity analysis for multivalued treatment
settings. We report the percentage average bias (in SD units) in the lower and
upper bounds of the point estimate interval, the non-coverage rate of the
confidence interval, the median interval of point estimates, and the median
confidence intervals calculated from the $1000$ simulated datasets.

The simulation results under the proposed framework of sensitivity analysis
are presented in Table table. When there is adequate overlap (Scenario-I), we
observe that the average bias in the SIPW point estimators of the pairwise
ATEs lies within $10\%$ of its standard deviation in almost all cases under
under the proposed framework. The percentile bootstrap confidence intervals
also satisfy the nominal $90\%$ coverage rate in the presence of adequate
overlap. However, we observe that the performance of the SIPW point estimators
and the percentile bootstrap confidence intervals deteriorates when there is a
lack of overlap (Scenario-II). These findings are similar to the ones observed
by Zhao et al. (2019) and Basit et al. (2023), and therefore, we recommend to
interpret the sensitivity analysis results under the proposed framework with
caution when there is a lack of overlap in the covariate distribution among
different treatment groups.

##  section Real Data Application

In this section, we applied our proposed sensitivity analysis framework for
multivalued treatments to estimate the causal effect of maternal education on
the fertility rate in Bangladesh. The dataset used in this analysis was
obtained from the latest round of the Bangladesh Multiple Indicator Cluster
Survey (MICS) 2019 (BBS & UNICEF, 2019). We considered the number of children
ever born to women in the reproductive age ($15-49$ years) as a measure of
fertility, which is the outcome of interest. The treatment variable was the
maternal education level that has four levels- ”pre-primary or none”,
”primary”, ”secondary”, and ”higher secondary or beyond”. The following
variables were considered as the observed confounders, which were assumed to
be associated with both treatment and outcome: place of residence (”urban” or
”rural”), district, religion (”muslim” or ”non-muslim”), household’s wealth
index, women’s age at marriage, and education level of the household head
(Tomal et al., 2022). After discarding missing values from the outcome,
treatment, and the observed confounders, we obtained complete data on $53,481$
women aged $15-49$ years to conduct the analysis.

Using the proposed sensitivity analysis framework, we conduct the sensitivity
analysis for the three pairwise ATEs between ”pre-primary or none”, ”primary”,
and ”secondary” vs. ”higher secondary or beyond” level of maternal education,
denoted by $\tau_{1,4}$, $\tau_{2,4}$, and $\tau_{3,4}$, respectively. We
considered
$\Gamma_{0}=\exp(\gamma_{0})=\\{1,1.25,1.5,1.75,2.00,2.25,2.50,2.75\\}$ to
assess the sensitivity of the causal effect estimates to the presence of
varying magnitudes of unmeasured confounding. As discussed in Remark
section.theorem, we did not need to specify the values of the sensitivity
parameter $\Gamma_{1}$ to conduct the sensitivity analysis. Since the
treatment variable is ordinal, we fitted an ordinal logistic regression model,
namely the continuation ratio regression model, to estimate the generalized
propensity scores (GPSs) using the chosen observed confounders. The estimated
GPSs ranged from $0.01$ to $0.98$ with a mean of $0.25$ across the four
treatment groups. In light of our findings in the simulation studies in
Section section, we recommend to interpret the results obtained from the SIPW
estimators with caution as some of the estimated GPSs are close to zero.

Table table: SIPW point estimate intervals and 90% percentile bootstrap confidence intervals for the pairwise ATEs under the RR framework for different values of the sensitivity parameter $\Gamma_{0}$. The pariwise ATEs of pre-primary or no education, primary education, and secondary education vs. higher seondary or beyond education are denoted by $\tau_{1,4}$, $\tau_{2,4}$, and $\tau_{3,4}$, respectively Estimand | $\Gamma_{0}$ | $\gamma_{0}$ |  Point estimate interval |  90% confidence Interval
---|---|---|---|---
| 1 | 0.00 | $(1.97,1.97)$ | $(1.92,2.02)$
| 1.25 | 0.22 | $(1.59,2.46)$ | $(1.52,2.53)$
| 1.5 | 0.41 | $(1.22,2.83)$ | $(1.14,2.90)$
| 1.75 | 0.56 | $(0.91,3.14)$ | $(0.83,3.21)$
| 2 | 0.69 | $(0.66,3.39)$ | $(0.59,3.46)$
| 2.25 | 0.81 | $(0.46,3.61)$ | $(0.38,3.68)$
| 2.5 | 0.92 | $(0.29,3.82)$ | $(0.22,3.90)$
$\tau_{1,4}$ | 2.75 | 1.01 | $(0.13,4.01)$ | $(0.03,4.09)$
| 1 | 0.00 | $(1.46,1.46)$ | $(1.43,1.48)$
| 1.25 | 0.22 | $(1.15,1.95)$ | $(1.10,1.99)$
| 1.5 | 0.41 | $(0.82,2.27)$ | $(0.77,2.32)$
| 1.75 | 0.56 | $(0.56,2.57)$ | $(0.50,2.63)$
| 2 | 0.69 | $(0.30,2.83)$ | $(0.24,2.88)$
| 2.25 | 0.81 | $(0.07,3.06)$ | $(0.02,3.11)$
| 2.5 | 0.92 | $(-0.12,3.24)$ | $(-0.17,3.28)$
$\tau_{2,4}$ | 2.75 | 1.01 | $(-0.28,3.39)$ | $(-0.33,3.44)$
| 1 | 0.00 | $(0.75,0.75)$ | $(0.72,0.80)$
| 1.25 | 0.22 | $(0.45,1.17)$ | $(0.39,1.23)$
| 1.5 | 0.41 | $(0.13,1.49)$ | $(0.08,1.54)$
| 1.75 | 0.56 | $(-0.13,1.75)$ | $(-0.18,1.80)$
| 2 | 0.69 | $(-0.33,1.97)$ | $(-0.38,2.01)$
| 2.25 | 0.81 | $(-0.49,2.13)$ | $(-0.54,2.17)$
| 2.5 | 0.92 | $(-0.62,2.27)$ | $(-0.67,2.31)$
$\tau_{3,4}$ | 2.75 | 1.01 | $(-0.73,2.38)$ | $(-0.78,2.41)$
Figure figure: Graphical representation of the sensitivity analysis results
for the three pairwise ATEs under the RR framework. The solid error bars are
the range of point estimates and the dashed error bars (together with the
solid bars) are the confidence intervals. The circles/triangles/squares are
the mid-points of the solid bars.

For each value of the sensitivity parameter $\Gamma_{0}$, we estimated the
partially identified point estimate intervals and obtained $100(1-\alpha)\%$
confidence intervals using percentile bootstrap with $B=1000$. Table table and
Figure figure represent the sensitivity analysis results under the risk-ratio-
based framework. We can observe that in the case of the pairwise ATE between
”pre-primary or none” and ”higher secondary or beyond” level of education
($\tau_{1,4}$), the confidence intervals do not contain the null value of
zero, suggesting a significant non-zero causal effect, for at least
$\Gamma_{0}=2.75$. That is, the estimated pairwise ATE $\tau_{1,4}$ is
significantly different from zero even in the presence of unmeasured
confounders that cause the true (unobserved) GPSs to be 2.75 times higher than
the observed GPSs. Similarly, the ATE between the ”primary” and ”higher
secondary or beyond” level of education ($\tau_{2,4}$) is statistically
significant for at least $\Gamma_{0}=2.25$. However, the ATE between
”secondary” and ”higher secondary or beyond” education ($\tau_{3,4}$) is
insignificant for values of $\Gamma_{0}$ as small as about $1.5$.

Based on the conducted sensitivity analysis under our proposed framework, we
can imply that the estimated pairwise ATE between ”pre-primary or none” vs
”higher secondary or beyond” education ($\tau_{1,4}$) and ”primary” vs ”higher
secondary or beyond” education ($\tau_{2,4}$) are less sensitive to the
presence of unmeasured confounding in that very strong unmeasured confounders
are needed to invalid the causal estimates of $\tau_{1,4}$ and $\tau_{2,4}$.
This can be perceived as a substantial evidence of significant causal effect
of maternal education on the fertility of women in Bangladesh.

##  section Conclusion

Sensitivity analysis is crucial for the accurate interpretation of causal
conclusions drawn from observational studies. In this paper, we proposed a
risk-ratio-based modified marginal sensitivity model by extending the odds-
ratio-based sensitivity model of Zhao et al. (2019). We further extended our
proposed sensitivity model to observational studies with multivalued
treatments. As the sensitivity parameter in our proposed framework quantifies
the degree of unmeasured confounding using risk ratios of the observed and
true (generalized) propensity scores, it is easier and more intuitive to
interpret the sensitivity analysis results under the proposed framework
compared to that under odds-ratio-based frameworks. We illustrated that
estimation of the inverse probability weighting (IPW) causal effect estimators
under our proposed framework can be computed efficiently using a percentile
bootstrap approach. We also conducted simulation studies that suggest that our
proposed sensitivity analysis framework performs well when there is an
adequate overlap among the treatment groups. However, the performance
deteriorates due to the instability of the inverse probability weights when
there is a lack of overlap. Finally, we demonstrated sensitivity analysis
under the proposed framework using an empirical study, where we have estimated
the causal effect of maternal education on the female fertility in Bangladesh.

There are a number of further potential research directions. We are currently
working on incorporating other smooth causal effect estimators, such as the
doubly-robust augmented IPW (AIPW) estimators (Robins et al., 1994) and
generalized overlap weighting (GOW) estimators (Li et al., 2019) into our
proposed framework. Furthermore, we intend to work on calibrating the the
sensitivity parameter in our framework to the observed confounders. Zhang &
Small (2020) have recently worked on such calibration of the sensitivity model
of Gastwirth et al. (1998) in matched observational studies.

## References

  * Basit et al. (2023) Basit, M. A., Latif, M. A. & Wahed, A. S. (2023).  Sensitivity analysis for causal effects in observational studies with multivalued treatments.  arXiv preprint arXiv:2308.15986 .
  * BBS & UNICEF (2019) BBS & UNICEF (2019).  Bangladesh multiple indicator cluster survey 2019, progotir pathey.
  * Bind & Rubin (2021) Bind, M.-A. C. & Rubin, D. (2021).  The importance of having a conceptual stage when reporting non-randomized studies.  Biostatistics & Epidemiology 5, 9–18.
  * Bonvini et al. (2022) Bonvini, M., Kennedy, E., Ventura, V. & Wasserman, L. (2022).  Sensitivity analysis for marginal structural models.  arXiv preprint arXiv:2210.04681 .
  * Carnegie et al. (2016) Carnegie, N. B., Harada, M. & Hill, J. L. (2016).  Assessing sensitivity to unmeasured confounding using a simulated potential confounder.  Journal of Research on Educational Effectiveness 9, 395–420.
  * Charnes & Cooper (1962) Charnes, A. & Cooper, W. W. (1962).  Programming with linear fractional functionals.  Naval Research Logistics Quarterly 9, 181–186.
  * Cornfield et al. (1959) Cornfield, J., Haenszel, W., Hammond, E. C., Lilienfeld, A. M., Shimkin, M. B. & Wynder, E. L. (1959).  Smoking and lung cancer: recent evidence and a discussion of some questions.  Journal of the National Cancer Institute 22, 173–203.
  * Dorn et al. (2021) Dorn, J., Guo, K. & Kallus, N. (2021).  Doubly-valid/doubly-sharp sensitivity analysis for causal inference with unmeasured confounding.  arXiv preprint arXiv:2112.11449 .
  * Gastwirth et al. (1998) Gastwirth, J. L., Krieger, A. M. & Rosenbaum, P. R. (1998).  Dual and simultaneous sensitivity analysis for matched pairs.  Biometrika 85, 907–920.
  * Imbens (2000) Imbens, G. W. (2000).  The role of the propensity score in estimating dose-response functions.  Biometrika 87, 706–710.
  * Kallus et al. (2019) Kallus, N., Mao, X. & Zhou, A. (2019).  Interval estimation of individual-level causal effects under unobserved confounding.  In The 22nd international conference on artificial intelligence and statistics. PMLR.
  * Kang et al. (2007) Kang, J. D., Schafer, J. L. et al. (2007).  Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data.  Statistical Science 22, 523–539.
  * Li et al. (2019) Li, F. et al. (2019).  Propensity score weighting for causal inference with multiple treatments.  The Annals of Applied Statistics 13, 2389–2415.
  * Robins (2002) Robins, J. M. (2002).  Comment on ‘Covariance adjustment in randomized experiments and observational studies’.  Statistical Science 17, 309–321.
  * Robins et al. (1994) Robins, J. M., Rotnitzky, A. & Zhao, L. P. (1994).  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American Statistical Association 89, 846–866.
  * Rosenbaum (1987) Rosenbaum, P. R. (1987).  Sensitivity analysis for certain permutation inferences in matched observational studies.  Biometrika 74, 13–26.
  * Rosenbaum (2002) Rosenbaum, P. R. (2002).  Observational studies.  Springer.
  * Rosenbaum et al. (2002) Rosenbaum, P. R. et al. (2002).  Covariance adjustment in randomized experiments and observational studies.  Statistical Science 17, 286–327.
  * Rubin (1978) Rubin, D. B. (1978).  Bayesian inference for causal effects: The role of randomization.  The Annals of Statistics 6, 34–58.
  * Tomal et al. (2022) Tomal, J. H., Khan, J. R. & Wahed, A. S. (2022).  Weighted bayesian poisson regression for the number of children ever born per woman in bangladesh.  Journal of Statistical Theory and Applications 21, 79–105.
  * VanderWeele & Ding (2017) VanderWeele, T. J. & Ding, P. (2017).  Sensitivity analysis in observational research: introducing the e-value.  Annals of internal medicine 167, 268–274.
  * Yadlowsky et al. (2022) Yadlowsky, S., Namkoong, H., Basu, S., Duchi, J. & Tian, L. (2022).  Bounds on the conditional and average treatment effect with unobserved confounding factors.  The Annals of Statistics 50, 2587–2615.
  * Zhang & Small (2020) Zhang, B. & Small, D. S. (2020).  A calibrated sensitivity analysis for matched observational studies with application to the effect of second-hand smoke exposure on blood lead levels in children.  Journal of the Royal Statistical Society Series C: Applied Statistics 69, 1285–1305.
  * Zhao et al. (2019) Zhao, Q., Small, D. S. & Bhattacharya, B. B. (2019).  Sensitivity analysis for inverse probability weighting estimators via the percentile bootstrap.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 81, 735–761.