#  Probing Inner-Edge of Dead Zones in Protoplanetary Disks with ALMA and Next
Generation Very Large Array

Takahiro Ueda National Astronomical Observatory of Japan, Osawa 2-21-1,
Mitaka, Tokyo 181-8588, Japan  Luca Ricci  Department of Physics and
Astronomy, California State University Northridge, 18111 Nordhoff Street,
Northridge, CA 91330, USA  Mario Flock Max-Planck Institute for Astronomy,
Königstuhl 17, 69117 Heidelberg, Germany  Zachary Castro  Department of
Physics and Astronomy, California State University Northridge, 18111 Nordhoff
Street, Northridge, CA 91330, USA

###### Abstract

The discovery of substructures in protoplanetary disks with ALMA has provided
us key insights on the formation of planets. However, observational
constraints on the formation of rocky planets have been still sparse,
especially because of the limited spatial resolution. The inner edge of so-
called dead zone is one of the preferential sites of rocky planet formation.
We investigate the capabilities of ALMA and ngVLA for observing a dust
concentration expected at the inner edge of the dead-zone around a Herbig
star. Herbig Ae/Be stars are useful laboratories for exploring the evolution
of rocky grains in protoplanetary disks because of their high luminosity which
pushes the dead-zone inner edge outward. We find that, thanks to its
unprecedented angular resolution and sensitivity, ngVLA can detect the dust
concentration at the dead-zone inner edge, with a reasonable integration time
of 10 hrs at $\lambda=3,7$ mm and 1 cm. The dust concentration is expected to
be optically thick at the ALMA wavelengths and cannot be spatially resolved
due to its limited resolution. On the other hand, the flux density from the
inner disk regions ($\sim$3–4 au) observed with current VLA is higher for
disks with a dust ring, and hence would be a useful indicator that help us
choose potential candidates of disks having a dust concentration at the inner
most region. With these observations we can characterize the process of dust
concentration in the innermost disk regions, where rocky planets can form.

††journal: ApJ††software: RADMC-3D (Dullemond et al., 2012)

##  1 Introduction

Understanding the formation process of rocky planets is of great importance in
planetary science. One preferential site of rocky planetesimal/planet
formation is the inner edge of the so-called dead zone (Kretke et al. 2009;
Tan et al. 2015; Ueda et al. 2019; Flock et al. 2019; Jankovic et al. 2021;
Ueda et al. 2021). The dead zone is the location where magneto-rotational
instability (MRI, Balbus & Hawley 1998) is suppressed because of poor gas
ionization (Gammie, 1996). The dead zone is likely to have an inner edge where
the gas temperature $T$ reaches $\sim$800–1000 K, above which thermal
ionization of the gas is effective enough to activate MRI (Gammie, 1996; Desch
& Turner, 2015). Across the dead-zone inner edge, the turbulent viscosity
arising from MRI steeply decreases from inside out, resulting in a local
maximum in the radial profile of the gas pressure (e.g., Dzyurkevich et al.
2010; Lyra & Mac Low 2012; Flock et al. 2016; Flock et al. 2017). The pressure
maximum traps solid particles (Whipple, 1972; Adachi et al., 1976) and the
local dust-to-gas mass ratio increases, leading potentially to the formation
of rocky planetesimals via the streaming instability (Youdin & Goodman, 2005;
Johansen et al., 2007; Carrera et al., 2015; Abod et al., 2019) or via the
gravitational instability (Coradini et al., 1981).

It is definitely important to observationally probe whether dust grains indeed
pileup at the dead-zone inner edge. Until now, most of the observational
constraints from the inner disk are derived from the near infrared (NIR)
emission. First observations of young Herbig Ae/Be stars found a NIR excess,
which could not be explained by simple radiative disk models (Hillenbrand et
al., 1992). Since then, the theoretical models of the inner rim advanced,
including the effect of the inner rim shape (Dullemond et al., 2001), dust
halo (Vinković et al., 2006), grain sizes (Kama et al., 2009) and
hydrodynamics (Flock et al., 2016). However, the near infrared emission allows
us only to probe the very hot (above 1000 K) and low density regions of the
uppermost disk atmosphere (Lazareff et al., 2017; Gravity Collaboration et
al., 2019). To probe the planet formation in the disk midplane we need high
angular resolution observations at longer wavelengths.

In the past years, observations with the Atacama Large
Millimeter/submillimeter Array (ALMA) have shown that ring and gap structures
are very common in protoplanetary disks (e.g., ALMA Partnership et al. 2015;
Andrews et al. 2018; Andrews 2020). However, most of observed rings are
located far from the central star (stellocentric radii $>$ 10–20 au) and are
thought to be composed of icy grains. In our work we will for the first time
model and present predictions for the dust concentration regions at the inner
disk close to the dead-zone edge. In our work we will explore the potential of
current and future sub-mm and radio interferometers, such as ALMA and the Next
Generation Very Large Array (ngVLA, Selina et al., 2018), to observe the dust
emission from these regions in protoplanetary disks.

##  2 Models and Methods

###  2.1 Dust disk

In this work, we use the surface density model obtained by Ueda et al. (2019)
which investigated the evolution of gas and dust in a disk around a Herbig
star. The central star has a mass of $2.5~{}M_{\odot}$, radius of
$2.5~{}R_{\odot}$ and effective temperature of 10000 K, resulting in a
luminosity of $56~{}L_{\odot}$. Our model parameters are similar to those of,
e.g., HD179218 and BF Ori (63.1 and 61.6$L_{\odot}$, respectively; Mendigutía
et al. 2011). In the simulations, the evolution of rocky grains in a gas disk
with a dead-zone is computed by taking into account the collisional
growth/fragmentation, turbulent diffusion and radial drift. The dust-size
evolution is computed with a single-size approach in which we compute the
evolution of grains with a representative size, i.e., maximum size, in each
radial grid cell (Sato et al., 2016). As the initial condition, the dust
surface density is set to be 100 times lower than the gas surface density
which is determined by the assumption of a steady state accretion with an
accretion rate of $10^{-8}~{}M_{\odot}~{}\rm{yr}^{-1}$. In this work we use
snapshots at $t=3\times 10^{5}$ yr of the simulations. We refer the readers to
Ueda et al. (2019) for the details of the simulations.

Figure 1:  Dust surface density and dust size profile around the dead-zone
inner edge.

We consider three models for the dust distribution named “No-pileup”,
“Fiducial” and “Extreme”, as shown in Figure 1. The No-pileup model
corresponds to the disk with no dead-zone inner edge. In this model, the
turbulence parameter $\alpha_{\rm turb}$ is set to be $10^{-3}$ and is uniform
in the entire region. Because of the uniform turbulence, the No-pileup model
has a smooth radial density profile with no dust concentration. The dust
surface density steeply decreases at 0.55 au where rocky grains evaporate into
the gas phase. The Fiducial model has $\alpha_{\rm turb}=10^{-2}$ and
$10^{-3}$ in the MRI-active and MRI-dead regions, respectively, with a steep
transition at the dead-zone inner edge located at $\sim$ 1 au. In this model,
due to the relatively strong turbulence, the dust size is regulated to $<1$ cm
by the turbulence-induced collisional fragmentation.

Figure 2:  Extinction optical depths in the vertical direction for a face-one
oriented disk models at ALMA and ngVLA wavelengths. The lower x-axis shows the
radial distance and the upper x-axis corresponds to the angular separation
with a distance of 140 pc.

Owing to their small size, dust particles are marginally coupled with the gas
and hence the turbulent diffusion prevents dust from accumulating efficiently.
The midplane dust-to-gas mass ratio is $\sim$ 3, which would trigger the
formation of rocky planetesimals via the streaming instability. The Extreme
model is an extreme case where the dust particles strongly pile up at the
dead-zone inner edge. This model has $\alpha_{\rm turb}=10^{-3}$ and $10^{-4}$
in the MRI-active and MRI-dead regions, respectively. In this model, due to
the weak turbulence, turbulent diffusion no longer prevents dust from trapping
at the dead-zone inner edge. In all models, the critical fragmentation
velocity of rocky grains is set to 10 ${\rm m~{}s^{-1}}$. The choice of the
fragmentation velocity is based on the theoretical models of collisions of
silicate aggregates composed of $\sim$ 0.1 ${\rm\mu m}$ monomers (e.g., Wada
et al. 2013). The exact value of the fragmentation velocity depends on the
composition, monomer size and surface energy, and has large uncertainty (e.g.,
Blum & Wurm 2008; Homma et al. 2019; Steinpilz et al. 2019). With these three
models we cover the most extreme cases, from having no pressure bump and no
dust concentration until having a pressure bump with very efficient dust
trapping.

###  2.2 Radiative transfer simulation

In order to obtain predictions for the maps of the dust continuum emission for
the models described above, we perform the three-dimensional radiative
transfer simulations with the Monte Carlo radiative transfer code RADMC-3D
(Dullemond et al., 2012). At each radial grid cell, the vertically integrated
grain size distribution is assumed to be a truncated power-law distribution
ranging from 0.1 ${\rm\mu m}$ to $a_{\rm max}$ with the power-law index of
3.5. The maximum grain radius $a_{\rm max}$ is determined by the simulation
and shown in Figure 1. In the vertical direction, we assume a mixing-settling
equilibrium for each dust size (Youdin & Lithwick, 2007); larger grains settle
more to the mid-plane. The grain size distribution is divided into 7 bins,
$<0.3~{}{\rm\mu m}$, 0.3–3 ${\rm\mu m}$, 3–30 ${\rm\mu m}$, 30–300 ${\rm\mu
m}$, 0.3–3 ${\rm mm}$, 3–30 ${\rm mm}$ and 3–30 ${\rm cm}$. We adopt a
representative size for each size-bin ($0.1~{}{\rm\mu m}$, 1 ${\rm\mu m}$, 10
${\rm\mu m}$, 100 ${\rm\mu m}$, 1 ${\rm mm}$, 10 ${\rm mm}$ and 10 ${\rm cm}$
for bin of $<0.3~{}{\rm\mu m}$, 0.3–3 ${\rm\mu m}$, 3–30 ${\rm\mu m}$, 30–300
${\rm\mu m}$, 0.3–3 ${\rm mm}$, 3–30 ${\rm mm}$ and 3–30 ${\rm cm}$,
respectively), and compute a spatial distribution and opacities of each size-
bin using the representative size. The opacities are computed using Mie theory
for each size bin assuming spherical compact grains with the material density
of 3.0 ${\rm g~{}cm^{-3}}$ and the optical constants of amorphous silicate
(Jaeger et al., 1994). We first perform the iterative thermal Monte Carlo
simulations and obtain the temperature structure of the disk, and then perform
the imaging simulation using the obtained temperature structure (see Ueda et
al. 2019).

Figure 2 shows vertical extinction optical depths of our disk models at ALMA
and ngVLA wavelengths. The optical depth of the No-pileup model has a smooth
radial profile and is close to unity at ALMA wavelengths and lower than unity
at ngVLA wavelengths. The Fiducial model has a peak in the optical depth
profile at $\sim$ 1 au. The optical depth at the peak reaches $\sim 100$ at
$\lambda\lesssim 3$ mm, while it is down to $\sim 3$ at ngVLA Band 4
($\lambda=1.0$ cm). In the Extreme model, the dust ring is quite optically
thick ($>$100) and extends to $\sim$ 2 au.

###  2.3 Simulations of ALMA and ngVLA observations

The synthetic images for the dust continuum emission obtained in the previous
section were used to simulate the results of observations with ALMA and ngVLA.
This is done using the CASA software package111https://casa.nrao.edu/,
following the method outlined in Ricci et al. (2018) (see also Harter et al.,
2020; Blanco et al., 2021). Synthetic observations with ALMA were produced at
wavelengths of 870 ${\rm\mu m}$ and 1.3 mm, whereas for the ngVLA we focused
on wavelengths of 3 and 7 mm, and 1 cm. In order to provide a fair comparison
between the imaging capabilities of these two facilities, the Declination sky
coordinate of the disk center was set to $+24~{}\deg$ (Taurus region) and to
$-24~{}\deg$ (Ophiuchus region) for the ngVLA and ALMA simulations,
respectively. At the latitudes of these observatories, these correspond to
maximum source elevations of about 80 and 89 degrees for the ngVLA and ALMA,
respectively.

The model synthetic images are Fourier transformed and sampled at the $(u,v)$
points of the ALMA and ngVLA observations using the SIMOBSERVE task in CASA.
For ALMA, we considered the antenna position file alma.out28.cfg available in
the CASA software package, which contains the longest baselines for this
array, and simulated aperture synthesis for an integration time of 10 hours.
For the ngVLA, we consider the REV-B array configuration, which contains 214
antennas of 18-meter diameter. The interferometric visibilities are manually
corrupted using the SETNOISE task. The rms noise values obtained on the maps
correspond to an integration time of 10 hours according to the ngVLA array
design specified in Selina et al. (2018).

The imaging of the interferometric visibilities is performed using the TCLEAN
task. For both the ALMA and ngVLA datasets we adopt a Briggs weighting scheme
with robust parameter $R=-1$ and multiscale deconvolver. These choices provide
a good trade-off between angular resolution and sensitivity for the emission
scales and brightness of the disk models considered in this work (e.g., Blanco
et al., 2021). For example, imaging with robust parameters $R>-1$ provides
better signal-to-noise ratios but at the expense of angular resolution and
some of the substructures presented in this work are not spatially resolved;
conversely, imaging with robust parameters $R<-1$ provides better angular
resolution, but the lower signal-to-noise for the surface brightness does not
allow us to detect many of the fainter disk substructures predicted by our
models. The distance to all our disk models was set to 140 pc, which is close
to the average distance of several nearby star forming regions.

##  3 Results

###  3.1 Synthetic images and radial intensity profile

Figure 3:  Synthetic images of the innermost region of disks observed with
ALMA and ngVLA. Each row corresponds to the No-pileup (top), Fiducial (middle)
and Extreme (bottom) model. The integration time is set to 10 hrs. The
observing beam-size is denoted with a white-filled ellipse at lower left in
each panel.

Figure 3 shows the synthetic images of the inner regions of the disk models
considered here, as observed with ALMA and ngVLA. The angular resolutions of
these synthetic observations are: $10.6\times 9.74$ mas (ALMA at 0.87 mm),
$15.4\times 14.3$ mas (ALMA at 1.3 mm), $2.04\times 1.14$ mas (ngVLA at 3 mm),
$4.75\times 2.66$ mas (ngVLA at 7 mm), $6.30\times 3.51$ mas (ngVLA at 1 cm).
The corresponding radial profiles of the azimuthally averaged intensity are
shown in Figure 4. The peak intensity, rms noise and peak signal-to-noise
ratio (SNR) in each map are listed in Table 1.

As a result of the limited angular resolution, the ALMA Band 6 and 7
observations show smooth intensity images with no clear signature of the dust
concentration at the dead-zone inner edge. In the No-pileup model, a bright
ring at $r\sim 0.5$ au, corresponding to the inner rim of the dust disk, is
observed at ngVLA Band 6 with a SNR $\approx~{}50$. In ngVLA Band 4 and 5
images, however, the dust rim is not observed clearly. This is because the
disk is optically thin so that the rim-structure is not detected with
sufficient SNR at these wavelengths even after 10 hrs of integration time.

In the Fiducial and Extreme models, on the other hand, the ngVLA maps a bright
ring at $\sim$ 1 au is clearly detected, which is caused by the dust-pileup at
the dead-zone inner edge. In the Fiducial model, the ring is detected with SNR
ratios of 45.4, 29.3 and 27.5 at ngVLA Band 6, 5 and 4, respectively. In the
Extreme model, the ring is detected with SNR ratios of 82.3, 61.7 and 70.7 at
ngVLA Band 6, 5 and 4, respectively. These high SNR ratios indicate that the
dust ring at the dead-zone inner edge can be detected even with a
significantly shorter integration time. If we set the integration time to 3
hrs, the SNR would decrease by a factor of $\sqrt{3/10}\sim 0.5$, which is
still enough to detect the ring. In the Extreme model, both rings at the rim
of the dust disk and at the dead-zone inner edge are detected at ngVLA Band 6,
while they are not resolved at ngVLA Band 4 and 5 due to the poorer
resolution.

Figure 4:  Azimuthally-averaged radial intensity profile of the No-pileup model (left), Fiducial model (center) and Extreme model (right) observed at ALMA and ngVLA wavelengths. The radial bins are sampled at a width of 0.1 au. The error bar denotes the standard deviation of the intensities divided by the number of beam segments at each radial bin . The intensity is multiplied by a square of the wavelength to compensate for the spectral dependence of the Black-body radiation in Rayleigh-jeans limit.  Table 1: Peak intensity, rms noise and SNR of the dust emission in the synthetic images. | $\lambda$ | Peak | Noise | SNR
---|---|---|---|---
| [mm] | [mJy/beam] | [mJy/beam] |
| 0.87 | 1.32 | $7.68\times 10^{-3}$ | 172
| 1.3 | $8.89\times 10^{-1}$ | $4.48\times 10^{-3}$ | 198
No-pileup | 3.0 | $3.75\times 10^{-3}$ | $7.45\times 10^{-5}$ | 50.3
| 7.0 | $2.74\times 10^{-4}$ | $3.40\times 10^{-5}$ | 8.06
| 10 | $1.32\times 10^{-4}$ | $2.06\times 10^{-5}$ | 6.41
| 0.87 | 1.10 | $8.01\times 10^{-3}$ | 137
| 1.3 | $6.94\times 10^{-1}$ | $4.48\times 10^{-3}$ | 155
Fiducial | 3.0 | $3.49\times 10^{-3}$ | $7.68\times 10^{-5}$ | 45.4
| 7.0 | $9.96\times 10^{-4}$ | $3.40\times 10^{-5}$ | 29.3
| 10 | $5.66\times 10^{-4}$ | $2.06\times 10^{-5}$ | 27.5
| 0.87 | 1.95 | $7.80\times 10^{-3}$ | 250
| 1.3 | 1.63 | $4.45\times 10^{-3}$ | 366
Extreme | 3.0 | $6.78\times 10^{-3}$ | $8.24\times 10^{-5}$ | 82.3
| 7.0 | $2.11\times 10^{-3}$ | $3.42\times 10^{-5}$ | 61.7
| 10 | $1.47\times 10^{-3}$ | $2.08\times 10^{-5}$ | 70.7

###  3.2 Flux densities from the inner disk regions

As shown in the previous section, future observations with the ngVLA can
provide evidence for the dust concentration expected by models of dust
evolution at the dead-zone inner edge. Here we investigate the possibility to
detect a signature of the dust concentration at the dead-zone inner edge with
current VLA resolution.

Table 2: Flux density coming from inside 3.5 au | No-pileup | Fiducial | Extreme
---|---|---|---
| [mJy] | [mJy] | [mJy]
870 ${\rm\mu m}$ | 11.0 | 11.2 | 9.11
1.3 ${\rm mm}$ | 3.37 | 3.54 | 4.04
3.0 ${\rm mm}$ | $5.11\times 10^{-1}$ | $5.54\times 10^{-1}$ | $8.76\times 10^{-1}$
7.0 ${\rm mm}$ | $8.13\times 10^{-3}$ | $1.77\times 10^{-2}$ | $4.77\times 10^{-2}$
1.0 ${\rm cm}$ | $1.67\times 10^{-3}$ | $6.29\times 10^{-3}$ | $2.06\times 10^{-2}$

Table 2 summarizes the flux density coming from the inner 3.5 au region
(diameter of 7 au). The spatial resolution of 7 au corresponds to an angular
resolution of $\sim 50$ mas with the distance of 140 pc, which is comparable
to the highest resolution of current VLA. At the wavelength of 870 ${\rm\mu
m}$, the flux density has no significant difference between the No-pileup and
Fiducial model, while the Extreme model predicts lower flux density. This is
because the intensity in the Extreme model steeply drops radially behind the
dead-zone inner edge (see Figure 4) because the disk temperature is low due to
the shadowing by the dust-pileup (Ueda et al. 2019, see Appendix A). At
wavelengths longer than 1.3 mm, the Extreme model predicts higher flux density
than the other models because of higher optical depth. The Fiducial model has
flux density comparable to that of the No-pileup model at $\lambda=0.87,1.3$
and 3.0 mm. Interestingly, on the other hand, at $\lambda=7.0$ mm and 1 cm,
the Fiducial model predicts significantly higher flux density than the No-
pileup model because of the dust-ring at the dead-zone inner edge. Therefore,
VLA observations at $\lambda=7.0$ mm and 1 cm with the highest resolution
could be a useful indicator for inferring potential candidates of disks having
dust-concentration at the dead-zone inner edge, even though it cannot resolve
the ring-structure spatially. Here we note that the sub-cm/cm emission from
disks is often affected by emission from ionized gas via a variety of
different mechanisms (e.g., Ubach et al., 2017; Ricci et al., 2021). This
contribution can be estimated using observations at longer wavelength ($>1$
cm) where the emission is dominated by the ionized gas (Carrasco-González et
al., 2016).

##  4 Dead-zone inner edge in known Herbig Ae/Be disks

In this section, we discuss possible targets for the observation of the dead-
zone inner edge. In passive disks, the location of the dead-zone inner edge
$r_{\rm DZIE}$ can be roughly evaluated as (Ueda et al. 2017, see also Flock
et al. 2016)

$\displaystyle r_{\rm DZIE}\sim
0.8\left(\frac{L_{*}}{56~{}L_{\odot}}\right)^{1/2}~{}{\rm au}.$ (1)

Since the dead-zone inner edge is characterized as the region where the gas
temperature reaches $\sim$ 1000 K, disks around more luminous stars are more
suitable for probing the dead-zone inner edge. Using Equation (1), the
expected angular separation between the central star and the dead-zone inner
edge $\theta_{\rm DZIE}$ is evaluated as

$\displaystyle\theta_{\rm DZIE}\sim
5.7\left(\frac{L_{*}}{56L_{\odot}}\right)^{1/2}\left(\frac{d}{\rm
140~{}pc}\right)^{-1}{\rm mas}$ (2)

Figure 5 compares the expected angular separation of the dead-zone inner edge
of known 218 Herbig disks with the angular resolution of ALMA and ngVLA
obtained in our synthetic images. The stellar luminosities and distances are
taken from Vioque et al. (2018).

Figure 5:  Expected angular separation of the dead-zone inner edge for 218
Herbig disks. The stellar parameters (luminosity and distance) are taken from
Vioque et al. (2018). The vertical solid lines show the angular resolutions of
ALMA and ngVLA obtained in our synthetic images; 10, 15, 1.5, 3.6 and 4.7 mas
for ALMA Band 7, ALMA Band 6, ngVLA Band 6, ngVLA Band 5 and ngVLA Band4,
respectively. The horizontal dotted line represents the luminosity of the star
in our model.

We clearly see that if the angular resolution is 3 mas, which is comparable to
the best resolution of ngVLA, the dead-zone inner edge can be spatially
resolved for many disks around known Herbig Ae/Be stars.

Figure 6:  ALMA archival data of the MWC 297 disk observed at Band 6. The
synthesized beam-size is denoted with a gray-filled ellipse at lower right.
The expected location of the dead-zone inner edge is denoted with a gray-open
ellipse.

Even though disks around luminous Herbig Ae/Be stars are expected to be less
massive than those around T-Tauri star due to strong irradiation (Fuente et
al. 2003; Alonso-Albi et al. 2009), there are some disks which are found to
extend far beyond the dead-zone inner edge. As an example, the circumstellar
disk around MWC 297 would be one of the best targets for probing the dead-zone
inner edge. MWC 297 is a young pre-main-sequence Herbig Be type star with a
spectral type of B1.5V, a distance of $\sim$ 375 pc and a luminosity of $\sim
39000~{}L_{\odot}$ (Vioque et al., 2018). Owing to its high luminosity, the
dead-zone inner edge is expected to be located at $\sim$ 21 au ($\sim 56$ mas
for $d=375$ pc). The angular separation between the central star and the dead-
zone inner edge can be resolved into $\sim$ 37 beam segments at ngVLA Band 6.
Figure 6 shows the ALMA archival data of the MWC 297 disk observed at ALMA
Band 6 (2018.1.00814.S, PI: L. Maud). Although the angular resolution is not
high enough to probe the dead-zone inner edge, the previous observation
confirm the dust emission extended to $\sim$ 0.5 arcsec (188 au at $d=375$
pc). The total flux at 218 GHz ($\lambda=1.38~{}{\rm\mu m}$) is 273 mJy,
significantly brighter than most T-Tauri disks (Andrews et al., 2013; Ansdell
et al., 2018). Therefore, the MWC 297 disk would be one of the best targets to
explore the dust distribution around the dead-zone inner edge.

It is worth noting that the MWC 297 disk is expected to have a inner disk
component which extends inside of the dead-zone inner edge (Guzmán-Díaz et
al., 2021). It is very unclear if disks around Herbig Ae/Be stars extend
inward to near the the dead-zone inner edge because only a small number of
Herbig Ae disks has been mapped with high spatial resolution. From the
Spectral Energy Distribution (SED) analysis, some disks around Herbig Ae/Be
stars are expected to have gaps or inner cavities (Maaskant et al., 2013).
Therefore, it is necessary to select potential targets based on their SEDs.

##  5 Summary

We investigated the capabilities of ALMA and ngVLA to detect and spatially
resolve the concentration of dust expected by models of dust evolution at the
inner edge of the dead-zone around a Herbig star with luminosity of
$56~{}L_{\odot}$. Herbig Ae/Be stars are useful laboratories for exploring the
evolution of rocky grains because of their high luminosity. We used models of
dust evolution around the dead-zone inner edge given by Ueda et al. (2019) and
produced predictions for the dust continuum emission for future observations
with ALMA and ngVLA. We found that, thanks to its unprecedented angular
resolution and sensitivity, ngVLA can detect and spatially resolve a dust
concentration at the dead-zone inner edge, with a reasonable integration time
of 10 hrs at $\lambda=3,7$ mm and 1 cm. The dust concentration is expected to
be optically thick at the ALMA wavelength and cannot be spatially resolved due
to its limited resolution. On the other hand, we can potentially infer the
dust concentration with current VLA by using the flux density emitted from the
innermost region with a diameter of 7 au which is comparable to the highest
resolution of VLA with a distance of 140 pc. Finally, we examined the expected
location of the dead-zone inner edge for 218 Herbig Ae/Be stars and found that
we can spatially resolve the dead-zone inner edge for tens of those disks with
ngVLA resolution. Particularly, the MWC 297 disk would be one of the best
targets because the separation between the central star and the dead-zone
inner edge can be resolved into $\sim$ 37 beam segments with ngVLA Band 6.
These observations will allow us to shed light on the process of dust
concentration in the inner disk regions, an important process for the
formation of rocky planets/planetesimals in young protoplanetary disks.

We thank an anonymous referee for useful comments. T.U. is supported by JSPS
KAKENHI Grant Numbers JP19J01929. Part of this work was supported by NAOJ
Overseas Visit Program for Young Researchers, NINS. M.F. acknowledge funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation program (grant agreement No. 757957). This paper
makes use of the following ALMA data: ADS/JAO.ALMA$\\#2018.1.00814.{\rm S}$.
ALMA is a partnership of ESO (representing its member states), NSF (USA) and
NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ.

##  Appendix A Mid-plane temperature profiles

Figure 7 shows the mid-plane temperature profile obtained from RADMC-3D
simulations. Inside $\sim$ 0.5 au, the temperature follows a temp profile
given under the optically thin regime as most of the dust component
evaporates. At $\sim$ 0.5 au, the temperature steeply decreases because the
disk becomes optically thick for the stellar irradiation. For the models with
a dust pileup, the dust-pileup cast blocks off the stellar light and cast a
shadow just behind it, results in a dip in the temperature (Ueda et al.,
2019). In our models, the temperature between 0.5 and 1 au is lower than that
obtained by the simulations that treat a detailed rim structure (e.g.,
Dullemond et al. 2001; Flock et al. 2016, 2017; Ueda et al. 2017). The disk
structure around the disk rim sensitively depends on the dust sublimation and
direct irradiation onto the disk rim. The detailed modelling of the disk rim
is not trivial because the dust sublimation process and rim structure depends
on each other, and will be provided in a future work.

Figure 7:  Mid-plane temperature profiles of our models.

## References

  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192
  * Adachi et al. (1976) Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
  * Alonso-Albi et al. (2009) Alonso-Albi, T., Fuente, A., Bachiller, R., et al. 2009, A&A, 497, 117
  * Andrews (2020) Andrews, S. M. 2020, ARA&A, 58, 483
  * Andrews et al. (2013) Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41
  * Ansdell et al. (2018) Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
  * Blanco et al. (2021) Blanco, D., Ricci, L., Flock, M., & Turner, N. 2021, arXiv e-prints, arXiv:2108.00907
  * Blum & Wurm (2008) Blum, J., & Wurm, G. 2008, ARA&A, 46, 21
  * Carrasco-González et al. (2016) Carrasco-González, C., Henning, T., Chandler, C. J., et al. 2016, ApJ, 821, L16
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Coradini et al. (1981) Coradini, A., Magni, G., & Federico, C. 1981, A&A, 98, 173
  * Desch & Turner (2015) Desch, S. J., & Turner, N. J. 2015, ApJ, 811, 156
  * Dullemond et al. (2001) Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, , , ascl:1202.015
  * Dzyurkevich et al. (2010) Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70
  * Flock et al. (2016) Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ApJ, 827, 144
  * Flock et al. (2017) —. 2017, ApJ, 835, 230
  * Flock et al. (2019) Flock, M., Turner, N. J., Mulders, G. D., et al. 2019, A&A, 630, A147
  * Fuente et al. (2003) Fuente, A., Rodríguez-Franco, A., Testi, L., et al. 2003, ApJ, 598, L39
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 462, 725
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Perraut, K., Labadie, L., et al. 2019, A&A, 632, A53
  * Guzmán-Díaz et al. (2021) Guzmán-Díaz, J., Mendigutía, I., Montesinos, B., et al. 2021, A&A, 650, A182
  * Harter et al. (2020) Harter, S. K., Ricci, L., Zhang, S., & Zhu, Z. 2020, ApJ, 905, 24
  * Hillenbrand et al. (1992) Hillenbrand, L. A., Strom, S. E., Vrba, F. J., & Keene, J. 1992, ApJ, 397, 613
  * Homma et al. (2019) Homma, K. A., Okuzumi, S., Nakamoto, T., & Ueda, Y. 2019, ApJ, 877, 128
  * Jaeger et al. (1994) Jaeger, C., Mutschke, H., Begemann, B., Dorschner, J., & Henning, T. 1994, A&A, 292, 641
  * Jankovic et al. (2021) Jankovic, M. R., Owen, J. E., Mohanty, S., & Tan, J. C. 2021, MNRAS, 504, 280
  * Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022
  * Kama et al. (2009) Kama, M., Min, M., & Dominik, C. 2009, A&A, 506, 1199
  * Kretke et al. (2009) Kretke, K. A., Lin, D. N. C., Garaud, P., & Turner, N. J. 2009, ApJ, 690, 407
  * Lazareff et al. (2017) Lazareff, B., Berger, J. P., Kluska, J., et al. 2017, A&A, 599, A85
  * Lyra & Mac Low (2012) Lyra, W., & Mac Low, M.-M. 2012, ApJ, 756, 62
  * Maaskant et al. (2013) Maaskant, K. M., Honda, M., Waters, L. B. F. M., et al. 2013, A&A, 555, A64
  * Mendigutía et al. (2011) Mendigutía, I., Calvet, N., Montesinos, B., et al. 2011, A&A, 535, A99
  * Ricci et al. (2021) Ricci, L., Harter, S. K., Ercolano, B., & Weber, M. 2021, ApJ, 913, 122
  * Ricci et al. (2018) Ricci, L., Liu, S.-F., Isella, A., & Li, H. 2018, ApJ, 853, 110
  * Sato et al. (2016) Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15
  * Selina et al. (2018) Selina, R. J., Murphy, E. J., McKinnon, M., et al. 2018, in Astronomical Society of the Pacific Conference Series, Vol. 517, Science with a Next Generation Very Large Array, ed. E. Murphy, 15
  * Steinpilz et al. (2019) Steinpilz, T., Teiser, J., & Wurm, G. 2019, ApJ, 874, 60
  * Tan et al. (2015) Tan, J. C., Chatterjee, S., Hu, X., Zhu, Z., & Mohanty, S. 2015, ArXiv e-prints, arXiv:1510.06703
  * Ubach et al. (2017) Ubach, C., Maddison, S. T., Wright, C. M., et al. 2017, MNRAS, 466, 4083
  * Ueda et al. (2019) Ueda, T., Flock, M., & Okuzumi, S. 2019, ApJ, 871, 10
  * Ueda et al. (2021) Ueda, T., Ogihara, M., Kokubo, E., & Okuzumi, S. 2021, ApJ, 921, L5
  * Ueda et al. (2017) Ueda, T., Okuzumi, S., & Flock, M. 2017, ApJ, 843, 49
  * Vinković et al. (2006) Vinković, D., Ivezić, Ž., Jurkić, T., & Elitzur, M. 2006, ApJ, 636, 348
  * Vioque et al. (2018) Vioque, M., Oudmaijer, R. D., Baines, D., Mendigutía, I., & Pérez-Martínez, R. 2018, A&A, 620, A128
  * Wada et al. (2013) Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588


# Combining Doubly Robust Methods and Machine Learning for Estimating Average
Treatment Effects for Observational Real-world Data

Xiaoqing Tan
Department of Biostatistics, University of Pittsburgh, Pennsylvania, USA
and
Shu Yang
Department of Statistics, North Carolina State University, North Carolina, USA
and
Wenyu Ye, Douglas E. Faries, Ilya Lipkovich and Zbigniew Kadziola
Eli Lilly and Company, Indiana, USA  Yang is partially supported by the NSF
grant DMS 1811245, NIA grant 1R01AG066883, and NIEHS grant 1R01ES031651.

###### Abstract

Observational cohort studies are increasingly being used for comparative
effectiveness research to assess the safety of therapeutics. Recently, various
doubly robust methods have been proposed for average treatment effect
estimation by combining the treatment model and the outcome model via
different vehicles, such as matching, weighting, and regression. The key
advantage of doubly robust estimators is that they require either the
treatment model or the outcome model to be correctly specified to obtain a
consistent estimator of average treatment effects, and therefore lead to a
more accurate and often more precise inference. However, little work has been
done to understand how doubly robust estimators differ due to their unique
strategies of using the treatment and outcome models and how machine learning
techniques can be combined to boost their performance. Here we examine
multiple popular doubly robust methods and compare their performance using
different treatment and outcome modeling via extensive simulations and a real-
world application. We found that incorporating machine learning with doubly
robust estimators such as the targeted maximum likelihood estimator gives the
best overall performance. Practical guidance on how to apply doubly robust
estimators is provided.

Keywords: Augmented inverse probability weighting, Double score matching,
Penalized spline of propensity methods for treatment comparison, SuperLearner

##  1 Introduction

Randomized control trials (RCTs) are considered to be the gold standard for
establishing the causal effects of interventions. They evaluate interventions
among comparable groups. However, sometime it would be impossible to conduct
RCTs due to limited resources or ethical issues. Observational studies, on the
other hand, examine effects in “real world” settings without manipulation. As
there is no intervention, some individuals with certain characteristics may
have a different probability of being exposed to treatment than others,
meaning that the covariate information between treatment groups may be highly
imbalanced. Therefore, it’s important to adjust for covariate imbalance issues
in observational studies.

There are two ways of adjustment for observational studies. The first kind is
based on the treatment model, also known as the propensity score (PS) model
where the PS is defined to be the probability of being treated given
covariates. The common inverse propensity treatment weighted estimator falls
into this category. The idea of weighting is to create a weighted pseudo-
population where treatments are “randomized”. Another kind is based on the
outcome model. This outcome imputation approach tries to impute the missing
potential outcomes based on outcome modeling. Estimators based on PS modeling
require the correct treatment model, and estimators based on outcome modeling
require the correct outcome model. In practice, it’s common to use linear
models for PS and outcome modeling. This, however, could be problematic
because linearity may be an inappropriate assumption for PS and outcomes when
the response surface is nonlinear. To account for potential nonlinearity, more
flexible models are needed to be considered.

Doubly robust estimators combine the above two adjustments in a fortuitous way
that the causal estimator can be consistent if either the outcome model or the
treatment model is correctly specified. Recently, various doubly robust
estimators of different kinds such as weighting, matching, and regression have
been proposed in the literature to estimate the average treatment effects
(ATEs) (Glynn & Quinn, 2010; van der Laan & Rubin, 2006; Yang & Zhang, 2022;
Zhou et al., 2019). Although they all use the PS and outcome mean models, they
combine them differently. Also, doubly robust estimators especially the ones
derived from the semiparametric efficiency theory are known to have the rate
double robustness property (Kang & Schafer, 2007; Chernozhukov et al., 2018)
in the sense that they remain root-$n$ consistent and asymptotically normal
when using machine learning approaches to estimating the nuisance functions.
The unaddressed question is how different doubly robust estimators perform
coupled with machine learning approaches. Also, little work has been done to
understand the challenges of covariates selection, overlapping of covariate
distribution, and treatment effect heterogeneity for doubly robust estimators.

In this paper, we review multiple popular doubly robust methods from the
categories of matching, weighting, or regression, and compare their
performance using different PS and outcome modeling via extensive simulations
as well as a real-world application. We found that incorporating machine
learning with doubly robust estimators such as the targeted maximum likelihood
estimator gives the best overall performance on estimating ATEs. The main
contribution of the paper is that we conduct a comprehensive evaluation of the
empirical performance of estimators. Also, practical guidance in applying
these estimators is provided.

The remaining paper is organized as follows. In Section 2, we discuss both
doubly robust estimators and singly robust estimators in detail. Section 3
presents the extensive comparative simulations and Section 4 reports on
performance of these estimators on a real-world application. Section 5
provides practical guidance and concludes the paper.

##  2 Methodology: Singly and Doubly Robust Estimators

###  2.1 Notation, Assumptions, and Estimand

All aforementioned methods are based on the potential outcomes framework
(Neyman, 1923; Rubin, 1974). Let $X_{i}$ be the set of observed covariates,
$A_{i}$ be the binary treatment indicator, and $Y_{i}$ be the observed outcome
for subject $i=1,2,\ldots,n$. Let $Y_{i}(a)$ be the potential outcome had
subject $i$ been given treatment assignment $a$, where $a=1$ is the treatment
and $a=0$ is the control. We assume subjects are independent. The causal
estimand of interest is the average treatment effect $\tau$, which is defined
as $\tau=E\\{Y(1)-Y(0)\\}$.

In reality, since only one of the potential outcomes is observed and another
is missing, sometimes the fundamental problem of causal inference (i.e.
estimating average treatment effect) is posited as a missing data problem. To
estimate causal parameter $\tau$ from data with non-randomized treatment
assignment, the following causal assumptions are needed.

######  Assumption 1 (Stable Unit Treatment Value)

The potential outcomes of any individual are unrelated to the treatment
assignment of other individuals and there is no multiple versions of the
treatment. Sometimes this is referred to as “consistency” in that the observed
outcome is equal to the potential outcome under the actually assigned
treatment, i.e., the observed outcome is given by $Y=Y(A)=AY(1)+(1-A)Y(0)$.

######  Assumption 2 (Conditional unconfoundedness or treatment ignorability)

Given covariates $X$, treatment assignment is independent of potential
outcomes, i.e. $\\{Y(0),Y(1)\\}\perp A|X$.

######  Assumption 3 (Overlap or common support or positivity)

All subjects are possible to receive either arm of treatment, i.e.,
$0<P(A=1|X=x)<1\quad\forall x\in X$.

Assumption 1 ensures there is no interference among subjects and there is no
multiple versions of the treatment. For observational studies, because the
exposure to treatment is not controlled, treatment may be related to the way a
subject might potentially respond. Assumption 2 states that it may be possible
to identify all pre-treatment covariates that are predictors of treatment or
outcome. If $X$ contains all confounders, among subjects who share the same
$X$ there is no association between $A$ and potential outcomes. A common
practice, in reality, is to collect a large number of possible confounders in
order to mitigate the violation of this assumption. However, including all
available covariates in the analysis may introduce bias and variance of the
causal effect estimator (Pearl, 2011; Myers et al., 2011; Brookhart et al.,
2006; Yang et al., 2020). Variable selection is hence an important procedure
when estimating ATE. Assumption 3 adds a restriction on the joint distribution
of treatment assignment and covariates. Overlap is an important issue in
estimating treatment effects from non-randomized trials. It describes the
extent to which the range of data is the same across the two treatment groups.
In fact, lack of overlap may affect all types of estimators. For matching
estimators, that means it is difficult to find good matches; for weighting
estimators, small overlap can result in extremely large weights; and for
regression estimators, they may heavily rely on extrapolation. When Assumption
3 is violated, a common practice is to trim the sample to restrict inference
to the one with sufficient overlap (Yang & Ding, 2018), or to coarsen PS (Zhou
et al., 2015).

Given the above stated assumptions, the ATE can be identified from observed
data by conditioning on available covariates

$\displaystyle\tau=E\\{Y(1)\\}-E\\{Y(0)\\}=E\\{E(Y|A=1,X)-E(Y|A=0,X)\\}$ (1)

where the outer expectation is with respect to distribution of $X$ over entire
population.

In this section, we first review two kinds of singly robust estimators based
on either outcome modeling or treatment modeling, respectively in the sense
that the consistency of the estimators relies on the correctness of the
underlying model. Then we review various doubly robust estimators that combine
outcome modeling and treatment modeling in estimating ATE. Augmented inverse
probability treatment weighting (AIPTW, Lunceford & Davidian, 2004; Glynn &
Quinn, 2010; Cao et al., 2009; Robins et al., 2000) belongs to a class of
weighting estimators. AIPTW is a combination of the basic inverse probability
weighting estimator and a weighted average of the outcome imputation
estimators. AIPTW improves the IPW estimator by fully utilizing information
about both treatment and outcome. Targeted maximum likelihood estimation
(TMLE, van der Laan & Rubin, 2006; van der Laan & Rose, 2011) is a regression
estimator based on maximum likelihood estimation and includes a “targeting”
step that optimizes the bias-variance tradeoff for the causal estimand. Double
score matching (DSM, Yang & Zhang, 2022; Zhang et al., 2021) belongs to the
class of matching estimators. DSM matches on both propensity score and
prognostic score. Penalized spline of propensity methods for treatment
comparison (PENCOMP, Zhang & Little, 2009; Zhou et al., 2019) is an example of
doubly robust regression estimator. PENCOMP estimates causal effects by
imputing missing potential outcomes with flexible spline models using multiple
imputations.

###  2.2 Estimators for ATE Based on Outcome Modeling

A traditional way to adjust for covariate imbalance in observational studies
is via the formulation of a regression model for the outcome $Y$ on $A$ and
$X$. That is, we can estimate the regression $E(Y|A,X)$ by modeling on the
observed data. Given the stated assumptions Section 2.1, the ATE can be
identified by equation (1).

An example of an regression imputation estimator can be obtained by fitting a
linear regression of $Y$ given $A$ and all $X$, which is given by

$\displaystyle E(Y|A,X)=\alpha_{0}+\alpha_{A}A+X^{T}\alpha_{X}.$ (2)

Suppose this is indeed the true regression model, by equation (1),
$\tau=\alpha_{0}+\alpha_{A}\cdot 1+X^{T}\alpha_{X}-(\alpha_{0}+\alpha_{A}\cdot
0+X^{T}\alpha_{X})=\alpha_{A}.$ The ATE can be obtained directly from the
coefficient for $A$, i.e. $\alpha_{A}$ in the regression model. If the true
regression is specified, this estimator is consistent of $\tau$, that is,
$\hat{\alpha}_{A}\overset{p}{\to}\tau$. Hence, this regression imputation
estimator is a singly robust estimator in the sense that it is consistent when
the outcome model is correctly specified.

However, model (2) assumes a constant treatment effect and could be severely
biased in the case of heterogeneous treatment effects. In practice, treatment
effects may vary across subjects. The regression imputation estimator is
usually obtained by modeling outcome separately within each treatment arm
rather than by using a single model (2). Note that the regression above can be
made more general as a general parametric model, since $Y$ can be of any type.
The missing potential outcomes are then imputed by their predictions from the
corresponding posited models. The regression imputation estimator is given by
the difference in the averages of potential outcomes. This can help to address
heterogeneity in treatment effects, however, the issue of model
misspecification still exists. Also, in case of a near violation of Assumption
3, the outcome model-based approaches rely on extrapolation.

###  2.3 Estimators for ATE Based on Treatment Modeling

Another class of ATE estimator for covariate adjustment relies on the
treatment model, or the propensity score model. The propensity score is
defined as the probability of treatment given covariates, i.e.,
$e(X)=E(A|X)=P(A=1|X).$

Under the assumptions defined in Section 2.1, given the propensity score, the
potential outcomes and treatment assignment are independent (Rosenbaum &
Rubin, 1983), i.e., $\\{Y(0),Y(1)\\}\perp A~{}|~{}e(X)$. Traditionally, the
estimation of propensity is by using a logistic regression where
$e(X,\beta)=\operatorname{logit}^{-1}\\{\exp(\beta_{0}+X^{T}\beta_{X})\\}$.

Consider the inverse of the propensity score as a weight for the outcome,
under the assumptions stated in Section 2.1 and the true propensity score

$\displaystyle E\Big{\\{}\frac{ZY}{e(X)}\Big{\\}}$
$\displaystyle=E\Big{\\{}\frac{ZY(1)}{e(X)}\Big{\\}}=E\Big{[}E\Big{\\{}\frac{ZY(1)}{e(X)}\Big{|}Y(1),X\Big{\\}}\Big{]}$
$\displaystyle=E\Big{\\{}\frac{Y(1)}{e(X)}E(Z|Y(1),X)\Big{\\}}=E\Big{\\{}\frac{Y(1)}{e(X)}E(Z|X)\Big{\\}}=E\Big{\\{}\frac{Y(1)}{e(X)}e(X)\Big{\\}}=E\\{Y(1)\\}.$

Similarly, $E\\{(1-Z)Y/e(X)\\}=E\\{Y(0)\\}$.

A well-known common estimator based on propensity score is Inverse Probability
Treatment Weighting (IPTW) (Lunceford & Davidian, 2004). Specifically, IPTW
estimates $\tau$ by the difference of inverse probability of treatment
weighted averages, which is given by

$\displaystyle\hat{\tau}_{IPTW}=\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}Y_{i}}{\hat{e}(X_{i})}-\frac{1}{n}\sum_{i=1}^{n}\frac{(1-A_{i})Y_{i}}{1-\hat{e}(X_{i})}.$

The inverse weights create a pseudo-population where there is no confounding
so the weighted averages can reflect averages in the target population. If the
true model for the propensity score model is specified, $\hat{\tau}_{IPTW}$ is
a consistent estimator of $\tau$. Hence, IPTW is a singly robust estimator in
the sense that it is consistent when the treatment model is correctly
specified. A major drawback of the IPTW estimator is that IPTW is highly
unstable due to weighting by the inverse of the propensity score. If any
propensity score is close to zero or one, the IPTW estimator may be extreme.

Variable selection is an important consideration when constructing the
propensity score. Including predictors of treatment but not outcome, i.e.,
instrument variables, in the treatment model or the outcome model may amplify
bias and variance of the causal estimator (Pearl, 2011; Myers et al., 2011).
Including outcome predictors, on the other hand, could boost efficiency
(Brookhart et al., 2006; Yang et al., 2020; Tan, Abberbock, Rastogi & Tang,
2022). Therefore, variable selection is needed before the estimation of
treatment effects to remove variables not related to outcomes. Besides,
weighting estimators are inferior in the case of extreme propensity scores
(Kang & Schafer, 2007). Poor overlap in propensity score distributions can
result in extremely large weights, leading to an unstable estimator with a
large variance. Furthermore, model misspecification of the propensity score
would also lead to a biased causal estimate.

###  2.4 Augmented Inverse Probability Treatment Weighted (AIPTW)

AIPTW estimator is a weighting based estimator that improves IPTW by fully
utilizing information about both the treatment assignment and the outcome
(Glynn & Quinn, 2010). It is a combination of IPTW estimator and a weighted
average of the outcome imputation estimators. Specifically, AIPTW is given by

$\displaystyle\hat{\tau}_{AIPTW}=\frac{1}{n}\sum_{i=1}^{n}\Bigg{\\{}\Big{\\{}\frac{A_{i}Y_{i}}{\hat{e}(X_{i})}-\frac{(1-A_{i})Y_{i}}{1-\hat{e}(X_{i})}\Big{\\}}-\frac{A_{i}-\hat{e}(X_{i})}{\hat{e}(X_{i})\\{1-\hat{e}(X_{i})\\}}\Big{[}\\{1-\hat{e}(X_{i})\\}\hat{m}_{1}(X_{i})+\hat{e}(X_{i})\hat{m}_{0}(X_{i})\Big{]}\Bigg{\\}}$
(3)

where $m_{1}(X)$ is a postulated model for $E(Y|A=1,X)$ and $m_{0}(X)$ is a
postulated model for $E(Y|A=0,X)$. The first line of equation (3) is the same
as $\hat{\tau}_{IPTW}$ and the rest adjusts this estimator by a weighted
average of the two outcome imputation estimators. Rearranging terms in
equation (3), $\hat{\tau}_{AIPTW}$ can be given by

$\displaystyle\hat{\tau}_{AIPTW}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\Big{\\{}\frac{A_{i}Y_{i}}{\hat{e}(X_{i})}-\frac{A_{i}-\hat{e}(X_{i})}{\hat{e}(X_{i})}\hat{m}_{1}(X_{i})\Big{\\}}-\frac{1}{n}\sum_{i=1}^{n}\Big{\\{}\frac{(1-A_{i})Y_{i}}{1-\hat{e}(X_{i})}+\frac{A_{i}-\hat{e}(X_{i})}{1-\hat{e}(X_{i})}\hat{m}_{0}(X_{i})\Big{\\}}$

AIPTW is a doubly robust estimator in that as long as either the outcome model
is correct or the propensity score model is correct, $\hat{\tau}_{AIPTW}$ is a
consistent estimator for $\tau$ (Glynn & Quinn, 2010). Also, it enjoys good
large-sample theoretical properties that it can be shown to be asymptotically
normally distributed via derivation through the theory of $M$-estimation.
Bootstrap can also be used to obtain the standard error estimates (Imbens,
2004). These standard error estimators tend to be reasonable unless the
estimated propensity scores are very extreme as weighting estimators are
inferior in the case of extreme propensity scores (Kang & Schafer, 2007). In
those scenarios, AIPTW is less robust to data sparsity and near violation of
the positivity assumption (Glynn & Quinn, 2010).

Recently, machine learning has gained popularity in the field of causal
inference (Peters et al., 2017; Prosperi et al., 2020; Tan, Chang, Zhou &
Tang, 2022; Tan, 2023; Tan, Qi, Seymour & Tang, 2022). Chernozhukov et al.
(2018) shows that the regression imputation and IPTW estimators using machine
learning nuisance function estimators tend to have large finite sample biases.
AIPTW, derived from the semiparametric efficiency theory, on the other hand,
enjoys the rate double robustness when combined with machine learning (Kang &
Schafer, 2007; Chernozhukov et al., 2018). That is, these doubly robust
estimators remain root-$n$ consistent and asymptotically normal when using
machine learning approaches to estimating the nuisance functions.
Incorporating AIPTW with machine learning could help to mitigate the impact of
regularization bias and overfitting on causal estimate (Chernozhukov et al.,
2018).

###  2.5 Targeted Maximum Likelihood Estimation (TMLE)

TMLE, introduced by van der Laan & Rubin (2006), is a maximum likelihood-based
estimator that incorporates a “targeting” step that optimizes the bias-
variance tradeoff for the targeted estimator, i.e., ATE. Specifically, TMLE
obtains initial outcome estimates via outcome modeling and propensity scores
via treatment modeling, respectively. These initial outcome estimates are then
updated to reduce the bias of confounding, which generates the so-call
“targeted” predicted outcome values.

The estimation steps of TMLE are given as follows. First, initial outcome
estimates are constructed by $\hat{Y}_{1}=\hat{E}(Y|A=1,X)$ and
$\hat{Y}_{0}=\hat{E}(Y|A=0,X)$, respectively, and the propensity scores
$\hat{e}(X)$ are estimated through treatment modeling. Then the targeted steps
begin by first calculating the inverse propensity $H_{a}$ for each subject,

$\displaystyle
H_{1}(A=1,X)=\\{\hat{e}(X)\\}^{-1}\text{~{}~{}and~{}~{}}H_{0}(A=0,X)=-\\{1-\hat{e}(X)\\}^{-1}.$

This is similar in form to inverse probability weights. Then, for the
treatment arm and the control arm, separately, the observed outcome $Y$ is
regressed on those estimated inverse propensity with fixed intercepts. Take a
binary outcome as an example. A logistic transform can be applied with a
binary outcome $Y$ where

$\displaystyle
logit\\{E(Y|A=1,X)\\}=logit(\hat{Y}_{1})+\epsilon_{1}H_{1}\text{~{}~{}and~{}~{}}logit\\{E(Y|A=0,X)\\}=logit(\hat{Y}_{0})+\epsilon_{0}H_{0}.$

In this way, we are able to generate updated, or so-called “targeted”
estimates of the set of potential outcomes, incorporating information from the
treatment mechanism in order to reduce the bias. The predicted outcomes are
then updated to be

$\displaystyle
logit(\hat{Y}_{1}^{*})=logit(\hat{Y}_{1})+\hat{\epsilon}_{1}H_{1}\text{~{}~{}and~{}~{}}logit(\hat{Y}_{0}^{*})=logit(\hat{Y}_{0})+\hat{\epsilon}_{0}H_{0}$

The final estimates is given by calculating ATE as mean difference in targeted
predicted outcome pairs across individuals

$\displaystyle\hat{\tau}_{TMLE}=\frac{1}{n}\sum_{i=1}^{n}(\hat{Y}_{1i}^{*}-\hat{Y}_{0i}^{*}).$
(4)

The variance of the TMLE estimator is obtained based on the efficient
influence curve (Porter et al., 2011; van der Laan & Rubin, 2006; Díaz & van
der Laan, 2011; van der Laan et al., 2007). In general, TMLE and AIPTW are
both efficient and have the minimum asymptotic variance under the large-sample
theory. However, it has been shown that under finite sample size or
challenging scenarios such as misspecified models, and nearly violated
positivity, TMLE may still provide causal estimates in the range of ATE since
$\hat{Y}_{a}^{*}$ in equation (4) are range-preserving while AIPTW may not
(van der Laan & Rose, 2011; Porter et al., 2011; Pirracchio et al., 2015). The
estimation of TMLE is usually coupled with SuperLearner (Díaz & van der Laan,
2011; van der Laan et al., 2007) for $\hat{Y}_{a}$ and $\hat{e}(X)$, which is
an ensemble of multiple statistical and machine learning models. It learns an
optimal weighted average of those models by giving higher weights to more
accurate models, and has been proven to have high accuracy (Díaz & van der
Laan, 2011; van der Laan et al., 2007). The performance of TMLE continues to
boost with the help of SuperLearner. Note that hybrid estimators have been
proposed to resemble TMLE and AIPTW that use coarsened propensity score
estimates instead of model-based ones (Zhou et al., 2015). They may have
better performance in case of severe model misspecification. However, the
choice of coarsening mechanism and the coarsening parameters may introduce
extra challenges (Zhou et al., 2015), and thus these hybrid estimators are not
included for comparison in the later simulation studies.

###  2.6 Double Score Matching (DSM)

Matching methods are powerful because they can be used to re-create a
randomized trial that is hidden under the observational study. DSM is a
matching-based estimator that uses the double balancing properties of
propensity score $e(X)$ and prognostic score $\Phi(X)$, the latter obtained
via outcome modeling, before the matching is conducted (Yang & Zhang, 2022).
The prognostic score is defined as a balancing score where
$\\{Y(0),Y(1)\\}\perp A|\Phi(X)$ (Hansen, 2008). The combination of $e(X)$ and
$\Phi(X)$, the double score, is also shown to be a balancing score (Antonelli
et al., 2018). That is, $\\{Y(0),Y(1)\\}\perp A|\\{e(X),\Phi(X)\\}$.

DSM estimator enjoys the double robustness property in that this result holds
even if only one score is correctly specified. For unit $i$, the potential
outcome under $A_{i}$ is the observed outcome $Y_{i}$. The potential outcome
under $1-A_{i}$ is not observed but can be imputed by the observed outcomes of
the nearest $M$ units with $1-A_{i}$. Denote the augmented score
$S=\\{e(X),\Phi(X)^{T}\\}^{T}$ as the matching variable, $J_{S,i}$ as the
index set for these matched subjects for subject $i$, and
$K_{S,i}=\sum_{j=1}^{n}I(i\in J_{S,i})$ as the number of times subject $i$ is
used as a match. The matched set is constructed with distance metric such as
Mahalanobis distance on $S$ that combines propensity score and prognostic
score. The initial DSM estimator of $\tau$ is

$\displaystyle\hat{\tau}_{DSM}^{(0)}=\frac{1}{n}\sum_{i=1}^{n}(2A_{i}-1)(1+M^{-1}K_{S,i})Y_{i}.$

A de-biasing DSM estimator $\hat{\tau}_{DSM}$ suggested by Yang & Zhang (2022)
further corrects bias by the method of sieves. Correctly the bias may help to
improve finite sample performance in practice although this bias is
asymptotically negligible (Yang & Zhang, 2022). A wild bootstrap procedure is
used to obtain the confidence interval based on Otsu & Rai (2017) for matching
estimators (Yang & Zhang, 2022).

Matching methods tend to be more stable tools when the propensity score is
extreme (Stuart, 2010). Matching estimators are robust to model
misspecifications if the misspecified model belongs to the class of covariate
scores (Waernbaum, 2012). DSM is robust against model misspecification of
either the propensity score model or the prognostic score model (Antonelli et
al., 2018; Yang & Zhang, 2022). Specifically, DSM provides multiple
protections to model misspecification by positing multiple candidate models
for both propensity scores and prognostic scores. This helps DSM to achieve a
near nominal coverage even under model misspecification (Yang & Zhang, 2022).
Furthermore, DSM can serve as a dimensional reduction tool in high-dimensional
confounding. However, adding too many covariates could result in potential
bias as matching estimators may not work well on high-dimensional covariates
(Abadie & Imbens, 2006). It is also pointed out that the number of posited
models and their functional forms affect the efficiency of DSM in a complex
way, resulting in an unstable performance if there is a large number of
working models (Yang & Zhang, 2022; Zhao & Yang, 2021). Hence, variable
selection in matching estimators is needed to help identify outcome predictors
for better efficiency and remove bias in estimating the ATE. Zhang et al.
(2021) investigate the performance of DSM under different strategies of
variable selection, using a caliper, and matching with or without replacement,
providing the best practice. Also, as the success of matching depends on the
functional forms of posited models, flexible machine learning methods can be
adopted in the modeling of propensity scores and prognostic scores before
matching.

###  2.7 Penalized Spline of Propensity Methods for Treatment Comparison
(PENCOMP)

PENCOMP is a type of regression- and multiple imputation-based approach (Zhou
et al., 2019). It builds on the method of the penalized spline of propensity
prediction previously used in missing data problems (Little & An, 2004; Zhang
& Little, 2009). Specifically, PENCOMP obtains propensity score $e(X)$ via
treatment modeling and uses spline-based regressions with propensity score
included for outcome modeling. Under the assumptions in Section 2.1, PENCOMP
has a double robustness property for estimating ATE.

PENCOMP uses Rubin’s rules for combining multiply imputed datasets. The method
first generates a bootstrap sample $b$ from the original data stratified on
the treatment group. The propensity score $e(X)$ is then estimated. Then the
potential outcomes for the treatments not assigned to subjects are predicted
with regression models that include splines on the logit of the propensity to
be assigned that treatment, plus other outcome predictors. For each treatment
group, each regression model is fitted separately. Specifically, the
regression model fitted for subjects with $A_{i}=a$ is given by

$\displaystyle
E\\{Y_{i}(a)|X_{i},A_{i}=a\\}=s\\{\hat{e}(X_{i})\\}+g\\{\hat{e}(X_{i}),X_{i}\\},~{}~{}i\in\\{i:A_{i}=a\\}$

where $s\\{\hat{e}(X)\\}$ is a penalized spline with pre-specified knots, and
$g\\{\hat{e}(X),X\\}$ is a parametric function of outcome predictors as well
as the estimated propensity. The missing potential outcome of a subject is
then imputed based on the predictive distribution of $E\\{Y(a)|X,A\\}$. The
bootstrap estimate $\hat{\tau}_{PENCOMP}^{(b)}$ is the difference in the
treatment means based on the observed outcome and the imputed values of $Y$.
The above procedure is repeated multiple ($B$) times. The final ATE estimate
$\hat{\tau}_{PENCOMP}=B^{-1}\sum_{b=1}^{B}\hat{\tau}_{PENCOMP}^{(b)}$. The
confidence interval of $\hat{\tau}_{PENCOMP}$ is generated from this
procedure.

According to Zhou et al. (2019), PENCOMP achieves comparable performance with
AIPTW in terms of bias, RMSE, and coverage under settings of low confounding
and correctly specified model with linear settings. Also, PENCOMP has some
advantages in nonlinear settings compared to AIPTW (Zhou et al., 2019).
However, in the case of model misspecification, there exists a severe
overcoverage with wider confidence interval for PENCOMP compared to AIPTW even
under a linear setting Zhou et al. (2019). It is also pointed out by Kang &
Schafer (2007) that if regression models are misspecified, doubly robust
methods could suffer from larger bias compared to singly robust methods.
Besides, the choice of the splines and knots can be challenging in practice.

##  3 Simulation Studies

In this section, we design Monte Carlo simulations to compare each doubly
robust method under different scenarios mimicking real-world data where there
is high nonlinearity in the relationship between covariates and treatment, and
the relationship between covariates and outcome. Specifically, we consider
settings with complex data generative models with multivariate covariates. We
also consider cases of different degrees of separation of the propensity score
distributions where the propensity scores may be close to zero or one. In
Table 1, we provide the open-source software or code that implements the
surveyed DR estimators. The code of our simulation studies can be found on
GitHub (https://github.com/ellenxtan/RealWorld-DoublyRobustML).

Table 1: Open-source software or code that implements the surveyed DR estimators. DR estimators | Open source code or package
---|---
AIPTW |  AIPW (Zhong et al., 2021)
TMLE |  tmle (Gruber & van der Laan, 2012)
DSM |  dsmatch (Yang et al., 2020)
PENCOMP |  PENCOMP (Zhou et al., 2019)

###  3.1 Data Generating Process

We design the data generating procedure following the work of Leacy & Stuart
(2014), but with the nonlinearity of treatment and outcome surfaces
considered. The sample size is set to be $n=2000$ throughout. The relationship
of variables generated in simulation is illustrated in Figure 1. First of all,
we generate $m=9$ independent and identically distributed standard normal
variables $X\in R^{9}\sim N(0,1)$. Let $W\in R^{9}$ be a nonlinear
transformation of $X$ where $W_{1},W_{2},W_{3}$ are confounders (i.e.,
predictors of both treatment assignment and outcome), $W_{4},W_{5},W_{6}$ are
treatment (only) predictors or instrumental variables, $W_{7},W_{8},W_{9}$ are
outcome (only) predictors. Specifically, we let $W_{1}=\exp(X_{1}/2),$
$W_{2}=\exp(X_{2}/3),$ $W_{3}=X_{3}^{2},$ $W_{4}=X_{4}^{2},$ $W_{5}=X_{5},$
$W_{6}=X_{6},$ $W_{7}=X_{7}+X_{8},$ $W_{8}=X_{7}^{2}+X_{8}^{2},$
$W_{9}=X_{9}^{3}$. All $W$ are standardized to have mean zero and variance 1.

Figure 1: Illustration of variables involved in simulation studies.

We generate the binary treatment indicator $A$ following Bernoulli$\\{e(X)\\}$
where $e(X)$ is the propensity score model. The propensity score model is
designed to be a linear combination of confounders and treatment predictors in
terms of $W$. To assess the property of different degrees on the separation of
propensity score distributions, we consider two settings for the PS overlap to
test the performance of estimators. A large overlap is an ideal case, which
means there is a reasonable amount of common support for the treated subjects
and the control subjects. The case of a small overlap, on the other hand,
indicates a majority of the treated subjects may fail to find any suitable
control in their neighborhoods. This could be common in scenarios such as rare
diseases in practice. Specifically, we design the large overlap case to be
$logit^{-1}\\{e(X)\\}=(-3-W_{1}+2W_{2}-3W_{3}+3W_{4}+2W_{5}+W_{6})/15$. We
design the case of small overlap to be
$logit^{-1}\\{e(X)\\}=(-8W_{1}+1.5W_{2}+0.5W_{3}-0.5W_{4}+2.5W_{5}-0.5W_{6})/5$.

We generate a continuous outcome, which is a linear combination of confounders
and outcome predictors in terms of W. Specifically, we let
$Y(0)=-2+1.5W_{1}-2W_{2}+1.5W_{3}+2.5W_{7}-W_{8}+W_{9}+\epsilon$ where
$\epsilon\sim N(0,1).$ We assume that only $X$’s are included in the candidate
set therefore nonlinearity in the outcome and treatment models induced by
$W$’s should be explicitly modeled by the analyst. This mimics that in
practice there exists potential nonlinearity between covariates and outcome
and this nonlinearity is needed to be considered when fitting models. We
consider both a homogeneous treatment effect setting and a heterogeneous
treatment effect setting. For the homogeneous setting, we consider a constant
treatment effect. We use $\tau=0$ in our simulation. Specifically,
$Y(1)=Y(0)+\tau$. For the heterogeneous setting, we allow the treatment effect
to vary by some covariates, and here $W_{1}$ and $W_{3}$. Specifically,
$Y(1)=Y(0)+\tau+5W_{1}+3W_{3}+2W_{1}W_{3}$. We try $\tau=0$. The Observed
outcome $Y$ is given by $Y=Y(1)A+Y(0)(1-A)$.

To summarize, there are four scenarios in total, which are

  1. 1)

homogeneous treatment effect with a large PS overlap (Homo TE + large overlap)

  2. 2)

homogeneous treatment effect with a small PS overlap (Homo TE + small overlap)

  3. 3)

heterogeneous treatment effect with a large PS overlap (Hetero TE + large
overlap)

  4. 4)

heterogeneous treatment effect with a small PS overlap (Hetero TE + small
overlap)

###  3.2 Evaluation Metrics

We consider the following metrics averaged over observations and simulation
replications.

  * •

Bias: $bias=R^{-1}\sum_{r=1}^{R}\\{\hat{\tau}^{(r)}-\tau\\}$

  * •

MSE: $mse=R^{-1}\sum_{r=1}^{R}\\{\hat{\tau}^{(r)}-\tau\\}^{2}$

  * •

Confidence interval coverage:
$coverage=R^{-1}\sum_{r=1}^{R}\mathbbm{1}\\{\tau\in(\hat{\tau}^{(r)}_{L,0.05},\hat{\tau}^{(r)}_{U,0.05})\\}$

  * •

Confidence interval width:
$width=R^{-1}\sum_{r=1}^{R}\\{\hat{\tau}^{(r)}_{U,0.05}-\hat{\tau}^{(r)}_{L,0.05}\\}$

  * •

Type I error:
$\alpha=1-R^{-1}\sum_{r=1}^{R}\mathbbm{1}\\{0\in(\hat{\tau}^{(r)}_{L,0.05},\hat{\tau}^{(r)}_{U,0.05})\\}$

  * •

Variance ratio:
$var.ratio=R^{-1}\sum_{r=1}^{R}[var_{m}\\{\hat{\tau}^{(r)}\\}]\\{var_{b}(\hat{\tau})\\}^{-1}$
where $var_{m}(\hat{\tau}^{(r)})$ is squared estimated standard error
$\hat{\tau}^{(r)}$ for the $r$-th replication and $var_{b}(\hat{\tau})$ is
variance of $\hat{\tau}$ over $R$ replications. Variance ratio measures the
ratio between mean variance of an estimator over $R$ replicates and variance
of estimates from the $R$ replicates. It evaluates the performance of model-
based standard errors of $\hat{\tau}$ by comparing them with simulation
variance reflecting the true variability of estimated $\tau$.

###  3.3 Analysis Steps and Compared Estimators

For analysis steps, we generate simulated data of $n=2000$ subjects. We apply
the Lasso (Tibshirani, 1996), a variable selection technique before the
estimation in order to remove variables that are not related to outcomes. The
outcome predictors are chosen by using all $X$’s as covariates and the
observed outcomes as response. Five-fold cross-validation is used to select
the best tuning parameter in Lasso, with the cross-validation deviance within
1 standard error of the minimum, as recommended by Zhang et al. (2021).
Outcome predictors are obtained for each treatment arm respectively. We then
use GLM, GAM, or SuperLearner separately to model PS and/or outcome before
estimating the final ATE. The candidate learners in the SuperLearner library
are: linear regression, stepwise regression, GAM, and Bayesian Additive
Regression Trees (Chipman et al., 2010). Doubly robust estimators for
comparison are AIPTW, TMLE, DSM, and PENCOMP. We also include for comparison
the singly robust estimator IPTW, and a regression imputation estimator,
denoted as IMP, which fits a twin outcome model on covariates separated by the
treatment group and the control group and then imputes the missing potential
outcomes with the posited models. Note that IPTW involves fitting the PS model
only. We generate $R=1000$ replications of simulated data sets. For AIPTW and
TMLE, bootstrap CIs are also computed. Each replication uses $B=500$ for the
bootstrap CI.

(a)

(b)

Figure 2:  Performance of various causal estimators under the four simulation
scenarios. (a) Box plots of the estimators. The red dotted lines indicate the
ground truth ATE. (b) Performance of the estimators in terms of CI coverage
and width, type I error, and variance ratio. The red dotted lines indicate the
ideal value or threshold of corresponding metrics.

Figure 3: Summary of the performance of various causal estimators under the
four simulation scenarios. Colors indicate the performance of the estimator
with green, yellow, and red meaning good, average, and poor performance,
respectively.

###  3.4 Simulation Results

Figure 2 shows the performance of various causal estimators under the four
simulation scenarios. A summary of the performance of the estimators is given
in Figure 3. In general, doubly robust estimators outperform singly robust
estimators. IMP tends to have robust performance by chance under homogeneous
settings, but may fail under heterogeneous settings with GLM and GAM because
of extrapolation and model misspecification. With SuperLearner, IMP greatly
reduces bias and variance. IPTW has improved performance with GAM in terms of
bias and MSE. However, with GLM and SuperLearner, IPTW may suffer from large
bias and MSE. IPTW may achieve a small bias under homogeneous treatment
effects but fails under settings of heterogeneous treatment effects. Compared
to overlap between propensity score distribution, treatment effect
heterogeneity has a larger effect on the performance of doubly robust
estimators. Using SuperLearner for treatment and outcome modeling, doubly
robust estimators achieve the smallest bias and MSE. Using GLM for treatment
and outcome modeling could suffer from huge bias because the relationship
between covariates and the outcome is nonlinear. Using GAM would help improve
the bias and MSE (bias in particular), but some nonlinearity may still be hard
to capture. Under the relatively easy homogeneous treatment effect setting,
doubly robust estimators including AIPTW, TMLE, and DSM achieve a nominal
confidence interval. Under challenging settings such as the heterogeneous
treatment effect, all doubly robust estimators suffer from overcoverage.
Overall, TMLE and AIPTW enjoy the most favorable performance with
SuperLearner. They have minimal bias and MSE, especially with SuperLearner.
Under the relatively simple homogeneous treatment effect setting, TMLE and
AIPTW achieve a nominal coverage, and with SuperLearner, these two estimators
even achieve the smallest CI width and control type I error under 5%. Under
the challenging settings, SuperLearner could inflate the variance ratio and
show an overcoverage issue. With GAM, TMLE and AIPTW may be able to achieve a
nominal confidence interval by using bootstrap (see in the setting of
heterogeneous treatment effects with a small overlap). However, in terms of
model misspecification, TMLE is more robust than AIPTW with GLM and GAM. TMLE
tends to have a smaller bias and MSE. DSM is more robust in terms of bias in
the case of model misspecification compared to other doubly robust estimators
(typically revealed in settings of heterogeneous treatment effects). DSM with
GAM may outperform DSM with SuperLearner in terms of bias and variance ratio
(see in the setting of heterogeneous treatment effects with a large overlap).
PENCOMP suffers from severe bias and MSE with GLM or GAM. PENCOMP has improved
performance with SuperLearner in terms of bias and MSE. However, PENCOMP
suffers from a severe overcoverage and a large type I error in nearly all
settings.

##  4 A Real-World Application

We apply different causal estimators to a real-world application, the
Reflections study (REFL), which is a study of real-world examination of
fibromyalgia for longitudinal evaluation of costs and treatments (Robinson et
al., 2012). We focus the analysis on opioid treatment arm (OPI cohort), and
non-narcotic opioid-like treatment arm (TRA cohort). There are 544 patients in
total. The outcome of interest is the change from baseline to LOCF in the
total score of the Fibromyalgia Impact Questionnaire (FIQ), which is a
continuous variable ranging from 0 to 80. There are 69 covariates in total, 24
of which are continuous variables and the other 45 are binary variables.
Earlier studies showed there is no difference in FIQ among the treatment
groups (Robinson et al., 2012; Yang et al., 2016). Here we apply the compared
causal estimators to estimate the causal effect of treatments for fibromyalgia
on the FIQ score. The analysis steps are similar to those in Section 3.3.
Figure 4 shows the performance of different estimators using different PS and
outcome modeling. The estimated treatment effect is about 0.03 with confidence
intervals of all doubly robust estimators including zero, which indicates
there is no evidence that there are treatment effects between the OPI cohort
and the TRA cohort on the FIQ score. With SuperLearner, the estimators achieve
the smallest standard error.

Figure 4: Performance of various causal estimators on the real-data
application. Different colors imply different causal estimators, x-axis
differentiate the PS and/or outcome models. The red dotted line indicates a
zero ATE.

(a)

(b)

Figure 5:  Performance of various causal estimators on the simulated REFL
study considering both a homogeneous treatment effect design (left, Zero ATE)
and a heterogeneous treatment effect design (right, Nonzero ATE),
respectively. (a) Box plots of the estimators. The red dotted lines indicate
the ground truth ATE. (b) Performance of the estimators in terms of CI
coverage and width, type I error, and variance ratio. The red dotted lines
indicate the ideal value or threshold of corresponding metrics.

###  4.1 Simulated REFL Study

Motivated by the real-world REFL study, we are interested in evaluating the
effects of a higher dimension of covariates and the effects of binary
variables on estimating treatment effects. The simulated REFL study is hence
designed to mimic the data distribution in the real REFL study. We generate
covariates X using Iman-Conover transformation (Iman & Conover, 1982) to
simulate correlated covariates from the real REFL data. The simulated REFL
data is of sample size 2000 (1000 for OPI cohort and 1000 for TRA cohort).

To simulate the treatment indicator $A$, we first generate the treatment model
by fitting an XGBoost (Chen & Guestrin, 2016) with cross-validation tuning
using the real REFL data. This is referred to as the “true” treatment model.
Specifically, the treatment indicator is used as the outcome and all patient
characteristics are included as covariates. We then generate the simulated
treatment assignment $A$ for the simulated REFL data following
Bernoulli($\hat{p}$) with $\hat{p}$ estimated from the “true” treatment model.
Similarly, the simulated outcomes $Y$ are obtained from the outcome model
which is fit by an XGBoost with cross-validation tuning. This is referred to
as the “true” outcome model.

Two simulation scenarios are designed with different “true” outcome models.
The first one is the zero treatment effect (Zero ATE) scenario, where the
“true” outcome model is based on fitting an XGBoost to outcome data with no
treatment indicator, then using the predicted value from that model, denoted
as $\hat{Y}$, to generate simulated outcomes by adding a Gaussian noise with
variation obtained from the cross-validation process. The resulting simulated
outcomes are considered “observed” outcomes, and $\hat{Y}$’s are considered
“truth”. The second one is the nonzero treatment effect (Nonzero ATE)
scenario, where the “true” outcome model is constructed by fitting an XGBoost
to outcome data with $A$ simulated from the “true” treatment model as a
covariate and other covariates included in the candidate set. The
corresponding predicted value $\hat{Y}$ from that fitted model is used to
simulate data. There are 1000 simulated datasets generated, each one conducts
a bootstrap 500 times.

Figure 5 illustrate the performance of various causal estimators on the
simulated REFL study considering Zero ATE design and Nonzero ATE design,
respectively. Similar to findings in the simulation study in Section 3, doubly
robust estimators, in general, outperform singly robust ATE estimators. Using
SuperLearner for treatment and outcome modeling, doubly robust estimators
achieve the smallest bias and MSE compared to using GLM, or GAM. Overall TMLE
and AIPTW share similar performance. With GAM or SuperLearner, TMLE and AIPTW
may be able to achieve a nominal confidence interval by using bootstrap
confidence intervals. DSM has a relatively larger bias compared to TMLESL,
TMLE, and AIPTW. DSM may achieve a nominal confidence interval with GAM.
PENCOMP suffers from severe bias and MSE with GLM or GAM. PENCOMP has improved
performance with SuperLearner in terms of bias and MSE. However, PENCOMP
suffers from a severe overcoverage and a large type I error in nearly all
settings.

##  5 Practical Recommendations and Discussion

We have reviewed multiple doubly robust estimators and conducted simulations
across a broad range of data scenarios. We vary causal inference test settings
by adjusting a variety of knobs in the simulations, which include nonlinearity
of treatment and outcome surfaces, degree of overlap between treatment
distributions as well as treatment effect heterogeneity. We make use of a
powerful machine learning technique SuperLearner to help improve ATE
estimation. Also, various doubly robust estimators are applied to a real-life
application of fibromyalgia as an example. In particular, we find that
incorporating machine learning with doubly robust estimators such as the TMLE
gives the best overall performance. Although in general TMLE and AIPTW are
both efficient and have the minimum asymptotic variance under the large-sample
theory, under finite sample sizes TMLE tends to be more robust to data
sparsity and near violations of positivity assumption because of its range-
preserving procedure for the predicted outcome estimates. Similar finding has
been shown in previous studies such as van der Laan & Rose (2011); Porter et
al. (2011); Luque-Fernandez et al. (2018); Bahamyirou et al. (2019). DSM is
robust to model misspecification as a matching estimator, but tends to have a
larger MSE compared to TMLE and AIPTW. The regression-based PENCOMP shows the
least ideal performance among all doubly robust estimators in case of model
misspecification and challenging scenarios, even when pairing with
SuperLearner. Further research is needed to demystify the performance of
PENCOMP found in our simulation studies.

Our paper helps to provide guidelines for practical use of doubly robust
estimators. Based on our extensive and realistic simulations, we recommend to
estimate the ATE in the following steps:

  * •

Perform variable selection to select outcome predictors.

  * •

Model the PS and the outcomes with SuperLearner, separated by the treatment
group and the control group.

  * •

Estimate the ATE by applying TMLE with the estimated propensity and outcome
estimates.

  * •

Use bootstrap for variance estimation of the ATE.

This work has multiple limitations that should be noted. First, throughout the
paper we only consider Lasso for variable selection to illustrate the
importance of the variable selection procedure. There might be better ways to
remove treatment predictors such as using machine learning algorithms like
random forest and neural networks. Soft variable selection strategies may also
be used where the variable selection is conducted without requiring any
modeling on the outcome, and thus provides robustness against misspecification
(Tang et al., 2021). Secondly, our work focuses only on doubly robust methods,
which in part due to their advantages in robustness over traditional methods.
However, these methods still require correct specification of at least one of
the models. Recent research proposed the use of model averaging across many
methods to improve the robustness of comparative analyses (Zagar et al.,
2022). Future work should compare the operating characteristics of doubly
robust approaches to model averaging, or perhaps simply incorporating multiple
doubly robust methods within the model averaging framework. In addition, an
extension of evaluating the use of doubly robust estimators on survival data
could be explored in the future.

In summary, this work has provided best practice guidance on the use of doubly
robust methods for comparative analysis based on real world data. The use of
machine learning for variable selection and model development, along with
estimation of treatment effects using TMLE, are found to help improved
operating characteristics of doubly robust methods.

## References

  * (1)
  * Abadie & Imbens (2006) Abadie, A. & Imbens, G. W. (2006), ‘Large sample properties of matching estimators for average treatment effects’, Econometrica 74(1), 235–267.
  * Antonelli et al. (2018) Antonelli, J., Cefalu, M., Palmer, N. & Agniel, D. (2018), ‘Doubly robust matching estimators for high dimensional confounding adjustment’, Biometrics 74(4), 1171–1179.
  * Bahamyirou et al. (2019) Bahamyirou, A., Blais, L., Forget, A. & Schnitzer, M. E. (2019), ‘Understanding and diagnosing the potential for bias when using machine learning methods with doubly robust causal estimators’, Statistical methods in medical research 28(6), 1637–1650.
  * Brookhart et al. (2006) Brookhart, M. A., Schneeweiss, S., Rothman, K. J., Glynn, R. J., Avorn, J. & Stürmer, T. (2006), ‘Variable selection for propensity score models’, American Journal of Epidemiology 163(12), 1149–1156.
  * Cao et al. (2009) Cao, W., Tsiatis, A. A. & Davidian, M. (2009), ‘Improving efficiency and robustness of the doubly robust estimator for a population mean with incomplete data’, Biometrika 96(3), 723–734.
  * Chen & Guestrin (2016) Chen, T. & Guestrin, C. (2016), XGBoost: A scalable tree boosting system, in ‘Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining’, KDD ’16, Association for Computing Machinery, Association for Computing Machinery, New York, NY, USA, p. 785–794.
https://doi.org/10.1145/2939672.2939785

  * Chernozhukov et al. (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. & Robins, J. (2018), ‘Double/debiased machine learning for treatment and structural parameters’, The Econometrics Journal 21(1), C1–C68.
https://doi.org/10.1111/ectj.12097

  * Chipman et al. (2010) Chipman, H. A., George, E. I. & McCulloch, R. E. (2010), ‘BART: Bayesian additive regression trees’, The Annals of Applied Statistics 4(1), 266–298.
  * Díaz & van der Laan (2011) Díaz, I. & van der Laan, M. J. (2011), ‘Super learner based conditional density estimation with application to marginal structural models’, The International Journal of Biostatistics 7(1), 1–20.
https://doi.org/10.2202/1557-4679.1356

  * Glynn & Quinn (2010) Glynn, A. N. & Quinn, K. M. (2010), ‘An introduction to the augmented inverse propensity weighted estimator’, Political Analysis 18(1), 36–56.
  * Gruber & van der Laan (2012) Gruber, S. & van der Laan, M. (2012), ‘tmle: An r package for targeted maximum likelihood estimation’, Journal of Statistical Software 51, 1–35.
  * Hansen (2008) Hansen, B. B. (2008), ‘The prognostic analogue of the propensity score’, Biometrika 95(2), 481–488.
  * Iman & Conover (1982) Iman, R. L. & Conover, W.-J. (1982), ‘A distribution-free approach to inducing rank correlation among input variables’, Communications in Statistics-Simulation and Computation 11(3), 311–334.
  * Imbens (2004) Imbens, G. W. (2004), ‘Nonparametric estimation of average treatment effects under exogeneity: a review’, Review of Economics and Statistics 86(1), 4–29.
  * Kang & Schafer (2007) Kang, J. D. Y. & Schafer, J. L. (2007), ‘Demystifying double robustness: a comparison of alternative strategies for estimating a population mean from incomplete data’, Statistical Science 22(4), 523 – 539.
https://doi.org/10.1214/07-STS227

  * Leacy & Stuart (2014) Leacy, F. P. & Stuart, E. A. (2014), ‘On the joint use of propensity and prognostic scores in estimation of the average treatment effect on the treated: a simulation study’, Statistics in Medicine 33(20), 3488–3508.
  * Little & An (2004) Little, R. & An, H. (2004), ‘Robust likelihood-based analysis of multivariate data with missing values’, Statistica Sinica 14(3), 949–968.
http://www.jstor.org/stable/24307424

  * Lunceford & Davidian (2004) Lunceford, J. K. & Davidian, M. (2004), ‘Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study’, Statistics in Medicine 23(19), 2937–2960.
  * Luque-Fernandez et al. (2018) Luque-Fernandez, M. A., Schomaker, M., Rachet, B. & Schnitzer, M. E. (2018), ‘Targeted maximum likelihood estimation for a binary treatment: A tutorial’, Statistics in medicine 37(16), 2530–2546.
  * Myers et al. (2011) Myers, J. A., Rassen, J. A., Gagne, J. J., Huybrechts, K. F., Schneeweiss, S., Rothman, K. J., Joffe, M. M. & Glynn, R. J. (2011), ‘Effects of adjusting for instrumental variables on bias and precision of effect estimates’, American Journal of Epidemiology 174(11), 1213–1222.
  * Neyman (1923) Neyman, J. (1923), ‘Sur les applications de la thar des probabilities aux experiences agaricales: Essay des principle. Excerpts reprinted (1990) in English’, Statistical Science 5, 463–472.
  * Otsu & Rai (2017) Otsu, T. & Rai, Y. (2017), ‘Bootstrap inference of matching estimators for average treatment effects’, Journal of the American Statistical Association 112(520), 1720–1732.
  * Pearl (2011) Pearl, J. (2011), ‘Invited commentary: understanding bias amplification’, American Journal of Epidemiology 174(11), 1223–1227.
  * Peters et al. (2017) Peters, J., Janzing, D. & Schölkopf, B. (2017), Elements of causal inference: foundations and learning algorithms, The MIT Press.
  * Pirracchio et al. (2015) Pirracchio, R., Petersen, M. L. & Van Der Laan, M. (2015), ‘Improving propensity score estimators’ robustness to model misspecification using super learner’, American Journal of Epidemiology 181(2), 108–119.
  * Porter et al. (2011) Porter, K. E., Gruber, S., Van Der Laan, M. J. & Sekhon, J. S. (2011), ‘The relative performance of targeted maximum likelihood estimators’, The International Journal of Biostatistics 7(1), 1–31.
  * Prosperi et al. (2020) Prosperi, M., Guo, Y., Sperrin, M., Koopman, J. S., Min, J. S., He, X., Rich, S., Wang, M., Buchan, I. E. & Bian, J. (2020), ‘Causal inference and counterfactual prediction in machine learning for actionable healthcare’, Nature Machine Intelligence 2(7), 369–375.
  * Robins et al. (2000) Robins, J. M., Hernan, M. A. & Brumback, B. (2000), ‘Marginal structural models and causal inference in epidemiology’, Epidemiology 11(5), 550–560.
  * Robinson et al. (2012) Robinson, R. L., Kroenke, K., Mease, P., Williams, D. A., Chen, Y., D’Souza, D., Wohlreich, M. & McCarberg, B. (2012), ‘Burden of illness and treatment patterns for patients with fibromyalgia’, Pain Medicine 13(10), 1366–1376.
  * Rosenbaum & Rubin (1983) Rosenbaum, P. R. & Rubin, D. B. (1983), ‘The central role of the propensity score in observational studies for causal effects’, Biometrika 70(1), 41–55.
  * Rubin (1974) Rubin, D. B. (1974), ‘Estimating causal effects of treatments in randomized and nonrandomized studies’, Journal of Educational Psychology 66(5), 688–701.
  * Stuart (2010) Stuart, E. A. (2010), ‘Matching methods for causal inference: a review and a look forward’, Statistical Science 25(1), 1–21.
  * Tan (2023) Tan, X. (2023), ‘Causal inference under data restrictions’, arXiv preprint arXiv:2301.08788 .
  * Tan, Abberbock, Rastogi & Tang (2022) Tan, X., Abberbock, J., Rastogi, P. & Tang, G. (2022), Identifying principal stratum causal effects conditional on a post-treatment intermediate response, in B. Schölkopf, C. Uhler & K. Zhang, eds, ‘Proceedings of the First Conference on Causal Learning and Reasoning’, Vol. 177 of Proceedings of Machine Learning Research, PMLR, pp. 734–753.
  * Tan, Chang, Zhou & Tang (2022) Tan, X., Chang, C.-C. H., Zhou, L. & Tang, L. (2022), A tree-based model averaging approach for personalized treatment effect estimation from heterogeneous data sources, in K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu & S. Sabato, eds, ‘Proceedings of the 39th International Conference on Machine Learning’, Vol. 162 of Proceedings of Machine Learning Research, PMLR, pp. 21013–21036.
https://proceedings.mlr.press/v162/tan22a.html

  * Tan, Qi, Seymour & Tang (2022) Tan, X., Qi, Z., Seymour, C. W. & Tang, L. (2022), ‘RISE: Robust individualized decision learning with sensitive variables’, arXiv preprint arXiv:2211.06569 .
  * Tang et al. (2021) Tang, D., Kong, D., Pan, W. & Wang, L. (2021), ‘Variable selection for doubly robust causal inference’, arXiv preprint arXiv:2007.14190 .
  * Tibshirani (1996) Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society: Series B (Methodological) 58(1), 267–288.
  * van der Laan et al. (2007) van der Laan, M. J., Polley, E. C. & Hubbard, A. E. (2007), ‘Super learner’, Statistical Applications in Genetics and Molecular Biology 6(1), 1–23.
  * van der Laan & Rubin (2006) van der Laan, M. J. & Rubin, D. (2006), ‘Targeted maximum likelihood learning’, The International Journal of Biostatistics 2(1), 1043–1083.
  * van der Laan & Rose (2011) van der Laan, M. & Rose, S. (2011), Targeted Learning: Causal Inference for Observational and Experimental Data, Springer Series in Statistics, Springer New York.
https://books.google.com/books?id=RGnSX5aCAgQC

  * Waernbaum (2012) Waernbaum, I. (2012), ‘Model misspecification and robustness in causal inference: comparing matching with doubly robust estimation’, Statistics in Medicine 31(15), 1572–1581.
  * Yang & Ding (2018) Yang, S. & Ding, P. (2018), ‘Asymptotic inference of causal effects with observational studies trimmed by the estimated propensity scores’, Biometrika 105, 487–493.
  * Yang et al. (2016) Yang, S., Imbens, G. W., Cui, Z., Faries, D. E. & Kadziola, Z. (2016), ‘Propensity score matching and subclassification in observational studies with multi-level treatments’, Biometrics 72(4), 1055–1065.
  * Yang et al. (2020) Yang, S., Kim, J. K. & Song, R. (2020), ‘Doubly robust inference when combining probability and non-probability samples with high dimensional data’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 82(2), 445–465.
  * Yang & Zhang (2022) Yang, S. & Zhang, Y. (2022), ‘Multiply robust matching estimators of average and quantile treatment effects’, Scandinavian Journal of Statistics pp. 1–31.
https://onlinelibrary.wiley.com/doi/abs/10.1111/sjos.12585

  * Zagar et al. (2022) Zagar, A., Kadziola, Z., Lipkovich, I., Madigan, D. & Faries, D. (2022), ‘Evaluating bias control strategies in observational studies using frequentist model averaging’, Journal of Biopharmaceutical Statistics pp. 1–30.
  * Zhang & Little (2009) Zhang, G. & Little, R. (2009), ‘Extensions of the penalized spline of propensity prediction method of imputation’, Biometrics 65(3), 911–918.
  * Zhang et al. (2021) Zhang, Y., Yang, S., Ye, W., Faries, D. E., Lipkovich, I. & Kadziola, Z. (2021), ‘Practical recommendations on double score matching for estimating causal effects’, Statistics in Medicine 15, 1–25.
  * Zhao & Yang (2021) Zhao, H. & Yang, S. (2021), ‘Outcome-adjusted balance measure for generalized propensity score model selection’, arXiv preprint arXiv:2107.12487 .
  * Zhong et al. (2021) Zhong, Y., Kennedy, E. H., Bodnar, L. M. & Naimi, A. I. (2021), ‘Aipw: An r package for augmented inverse probability weighted estimation of average causal effects’, American Journal of Epidemiology .
https://doi.org/10.1093/aje/kwab207

  * Zhou et al. (2015) Zhou, J., Zhang, Z., Li, Z. & Zhang, J. (2015), ‘Coarsened propensity scores and hybrid estimators for missing data and causal inference’, International Statistical Review 83(3), 449–471.
  * Zhou et al. (2019) Zhou, T., Elliott, M. R. & Little, R. J. (2019), ‘Penalized spline of propensity methods for treatment comparison’, Journal of the American Statistical Association 114(525), 1–19.