Vol.0 (20xx) No.0, 000–000

11institutetext: CAS Key Laboratory of Optical Astronomy, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012,
China<EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: Technical Institute of Physics and Chemistry, Chinese Academy
of Sciences, Beijing 100190, China
33institutetext: University of Chinese Academy of Sciences, Beijing 100049,
China
44institutetext: Department of Physics, University of California, Santa
Barbara, CA 93106, USA
Received 20xx month day; accepted 20xx month day

# Numerical study of a 20W class QCW pulsed sodium guide star laser’s
performances at five sites in China

Hong-Yang Li 1133 Lu Feng 11 Jun-Wei Zuo 22 Qi Bian 22 Bo-Tian Sun 44 Sui-Jian
Xue 11 Zhi-Xia Shen 11 Yang-Peng Li 11 Yong Bo 22

###### Abstract

In the past few years, Chinese astronomical community is actively testing
astronomical sites for several new optical/infrared ground-based telescopes.
These site testing campaigns conducted were mainly focused on fundamental
performances of the site, such as cloud coverage, seeing, temperature, etc.
With increasing interests in sodium laser guide star adaptive optics for these
new telescopes in the Chinese astronomical community, it is interesting to
investigate the performance of the laser guide star at these sites, especially
considering that the sodium laser guide star’s on-sky performance is
significantly influenced by sites’ local performances, such as geomagnetic
field, sodium layer dynamics, density of air molecule, etc. In this paper, we
studied sodium laser guide star’s performance of a 20W class Quasi-CW pulsed
laser developed by TIPC with numerical simulation for five selected sites in
China.

###### keywords:

instrumentation: adaptive optics — methods: numerical — site testing

##  1 Introduction

As a subsystem of Adaptive Optics (AO) system, sodium laser guide star system
can provide a stable and reliable reference source for telescopes’ AO systems
to correct atmospheric turbulence (Hardy 1998). Sodium laser guide star
systems have been utilized in most of those highly-productive large ground-
based telescopes , such as the Keck Telescope (Chin et al. 2016), the European
Southern Observatory’s (ESO) Very Large Telescope (Calia et al. 2014), the
twin Gemini Telescopes, the Gran Telescopio CANARIAS (GTC) et cetera, and such
system is also designed in the first place for those next generation 30m class
ground-based telescopes, such as the Thirty Meter Telescope, the European
Extremely Large Telescope (E-ELT) and the Giant Magellan Telescope (GMT).

One of the most important future astronomical projects within the Chinese
astronomical community, the 12m diameter Large Optical/infrared Telescope
(LOT) is also planned to have AO system installed as soon as the telescope is
constructed (Feng et al. 2020). By far, several excellent candidate sites for
this large telescope and several smaller telescopes to be built in near future
have been investigated in the western part of China, such as Tibet, Xinjiang,
Qinghai, Sichuan, etc (Feng & Hao 2020). Site monitoring scientists have made
great efforts in evaluating fundamental astronomical performances of these
sites, including cloud coverage, seeing, meteorological parameters,
preciptical water vapor, et cetera (Liu et al. 2012, 2015, 2016; Yao et al.
2013; Qian et al. 2015; Wang et al. 2015; Wu et al. 2016), but performances
relating to sodium laser guide star’s on-sky performances at these candidate
sites, like the geomagnetic field, sodium layer dynamics and so on, are not
being monitored due to the difficulties and costs for such measurements.
Therefore, we investigated the performance of sodium guide star at certain
typical sites in China in advance using our model combining with data from
remote sensing satellites as well as other atmospheric models, to quantify the
influence of some parameters related to sites on the performance of sodium
guide star, such as geomagnetic field and sodium layer dynamics. The results
of this study may shed lights on the future application of sodium guide stars
in the western part of China.

For comparative simulation study for the on-sky performance of sodium LGS, we
selected five sites for this study: Ali, Lenghu, Muztagh Ata, Gaomeigu in
western China and Xinglong (an important site in eastern China) as a
comparison. Figure 1 shows the location of the five sites. The Ali site is
located in the southwest of Tibet autonomous region with an altitude of 5100m.
It is one of the three candidate sites for the 12-meter Large Optical/infrared
Telescope (LOT) (Yao et al. 2013). The Muztagh Ata site is also one candidate
site for LOT, which is located in the Pamir Plateau, in the southwest of
Xinjiang Province of China, and more information about the site can be seen
from Xu et al. (2020). The Lenghu site is the site for the future 2.5m Wide
Field sky Survey Telescope(WFST) and it is located in the northwest border of
Qaidam Basin with an average altitude of 2800m. The Gaomeigu site is located
in Lijiang, Yunnan province with an altitude of 3200m, which has the largest
general-purpose optical telescope, Lijiang 2.4-meter Telescope (LJT), in
China. Xinglong site, located at the highest latitude among the selected
sites, has a history of more than 50 years and possesses two important
astronomical observation instruments, the LAMOST and the 2.16m Telescope.

The laser that we chose to simulate in this study is a 20W Quasi Continuous
Wave (QCW) pulsed laser from the Technical Institute of Physics and Chemistry
(TIPC) with whom we have been cooperating both in experiments and simulation
study of sodium laser guide star for years. Significant amount of lab/field
tests data have been gathered (Wei et al. 2012; Jin et al. 2014; Otarola et
al. 2016; Bian et al. 2016a, b, 2017, 2020), and were used to verify our model
previously (Jin et al. 2015; Feng et al. 2016). Although we are confident that
our model would be suitable to predict LGS generation performances for other
types of sodium laser, such as Continuous Wave (CW) laser, Macro-micro pulsed
laser (Li et al. 2021), we would like to present prediction results for such
lasers at different sites until sufficient field tests results have been
gathered to verify our model.

In the rest of the paper, we will describe parameters related to the laser and
sites in Section 2, and the simulation results will be presented in Section 3.

Figure 1: The location of the five sites in China (Xinglong, Lenghu, Muztagh
Ata, Ali and Gaomeigu)

##  2 METHOD

In this paper, we used a simulation tool called Photon Return Simulation
program (PRS) (Feng et al. 2015; Jin et al. 2015; Li et al. 2021) to predict
the on-sky performance of 20W Quasi-CW pulse laser at aforementioned sites in
China. The simulation is a Monte Carlo simulation based on rate equation.
During the simulated time period (for example 1s), the state of each atom
under the influence of laser light field as well as other different physics
such as geomagnetic field influence, light pressure, D2b repumping are being
tracked for every 30 nanosecond. The sampling time is short enough to simulate
a 100$\mu$s QCW laser pulse. We simulate this one second process repeatedly
for about 100 thousand times, which is sufficient for the result of the Monte
Carlo simulation to converge and to calculate the mean value of coupling
efficiency as in the Lidar equation as a representative for each laser power
density. After that, we could obtain the coupling efficiency versus the
projected on-sky power density curve. In this section, we will introduce the
sources and values of these parameters.

###  2.1 The parameters of TIPC laser

The spectrum and pulse format of sodium laser plays an important role in
predicting the return flux of sodium guide star (Kibblewhite 2000). As
described previously, the laser investigated in our study is a 20W quasi-CW
pulse laser from TIPC. Main parameters of the laser are summarized in Table 1.
The spectrum and single pulse format of the laser are shown in Figure 2. The
central wavelength of the laser is at sodium D2a line, 589.159nm, which
corresponds to 0 MHz in Figure 2. The Full Width Half Maximum (FWHM) of the
spectrum is about 0.3 GHz. The spectrum of the laser contains several
longitudinal modes and the distances between adjacent longitudinal modes are
equal which is about 225 MHz.

Parameter name | Value
---|---
Center wavelength |  $589.159$ nm
Linewidth |  $\sim 0.3$ GHz
Longitudinal mode interval |  $\sim 225$ MHz
Longitudinal mode width (FWHM) |  $5$ MHz
Pulse repetition rate |  $600$ Hz
Pulse width |  $\sim 100$ $\mu$s
Average power |  $20$ W
Polarization | circular
D2b repumping factor | $15\%$
Table 1: Laser parameters for TIPC 20 Watt Quasi-CW pulsed sodium laser

Figure 2: Spectrum format and single pulse format of TIPC 20W Quasi-CW pulse
laser. (left) The width (FWHM) of spectrum is 300MHz and the width (FWHM) for
each longitudinal modes is 5MHz. (right) The 100 $\mu s$ single pulse format
of the laser

###  2.2 Sodium atom column density

Photon return number is closely related to the sodium atom column density in
mesosphere, while the sodium Vertical Column Density(VCD) will change
obviously with latitude (Langowski et al. 2017). We utilize the empirical
formula of sodium column density proposed by Fussen, which is based on remote
sensing data from GOMOS (Fussen et al. 2004, 2010). The satellite is designed
for measuring column densities of molecules or atoms including sodium atom in
the atmosphere. Its results are further explored to show the variation of the
sodium VCD with different latitudes and different months. Based on these data,
Fussen formed an empirical formula of sodium column density against location
and month as shown in Equation 1. In this equation, N represents sodium
density for a certain month m (January is m=0) and latitude $\phi$ (in radians
for Eq. 1). Other parameters in Equation 1 are listed in Table 2. According to
this formula, we calculated the sodium column density in the mesosphere in
different months in China for the latitude ranging from 15°N to 55°N , and the
detailed results are visualized from Figure 3 (left). In general, the sodium
column density is higher in winter and lower in summer at the same latitude.
With the increase of latitude, it increases in winter and decrease in summer,
so the difference between the two seasons raises with latitude.

The sodium VCD reaches the highest in November, the lowest in May or June at
the selected sites in China, as shown in Figure 3 (right). The differences of
sodium VCD between winter and summer are relatively large. For example, even
for Gaomeigu site-the site with smallest annual difference among the sites-the
sodium VCD is 87.9% higher in November than in June: in November, it can reach
$4.1\times 10^{13}$ (atoms/$m^{2}$), whereas in June, it is only $2.2\times
10^{13}$ (atoms/$m^{2}$). Moreover, the difference of the sodium VCD between
different sites is relatively large from November to January next year, while
the difference between June and October is very small.

$\displaystyle
N(m,\phi)[cm^{-2}]=t_{0}+t_{1}\cos(\frac{2\pi}{12}m+t_{2})+t_{3}(\phi+\frac{\pi}{2})(\phi-\frac{\pi}{2})\cos(\frac{2\pi}{6}m+t_{4})$
$\displaystyle t_{i\leq 1}=f_{i}(a_{0}+a_{1}\phi+a_{2}\phi^{2}+a_{3}\phi^{3})$
$\displaystyle t_{0}=3.28\times 10^{9}$ (1)

i | $f_{i}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$
---|---|---|---|---|---
1 | $1\times 10^{9}$ | 0.1282 | 1.549 | 0.1780 | 0.03511
2 | 1 | 0.4017 | 0.8216 | -0.1282 | -0.2980
3 | $1\times 10^{9}$ | -0.2630 | 0.1121 | 0.6355 | -0.3566
4 | 1 | -1.5635 | -3.0526 | 1.3802 | 1.7637
Table 2: Parameters for Equation 1.

Figure 3: (left) Sodium column density changes with month and latitude (15
°N-55 °N) in China and is basically higher in winter and lower in summer.
(right) Five typical sites’ sodium column density varies with month. The
densities at Lenghu and Muztagh Ata are almost the same, so the curves almost
coincide.

###  2.3 The effect of geomagnetic field

The Larmor precession caused by the geomagnetic field will result in the M
state of sodium atoms to change periodically with time (Moussaoui et al. 2008,
2009). The period depends on the strength of the geomagnetic field and the
angle between the pointing direction of the laser and the geomagnetic field
vector. Since the transition cross section of the sodium atom is different
when it is at different M ground state, the periodical shift due to Larmor
precession will cause a degradation of the overall coupling efficiency for LGS
generation. The strength and direction of the geomagnetic field vector will
change with the location of the site. The coupling efficiency and photon flux
will change accordingly. This will lead to some differences in the coupling
efficiency and photon return number at these sites, so it is necessary to
determine the strength and direction of the geomagnetic field at these sites.
In Table 3, we listed the sites’ local geomagnetic field information,
including declination/inclination angle and the strength of the geomagnetic
field vertor, which is from the World Magnetic Model (WMM).

The pointing direction of laser will change as the telescope tracks the
observation target. Therefore, the angle between the laser and the geomagnetic
field will also change accordingly. In our study, we evaluate the on-sky
performance variation due to geomagnetic field effects. When we were analyzing
the on-sky performance with different projected laser power for different
months, we only simulate in the zenith direction because of two reasons: 1. to
easily compare with LGS field tests, since for these field tests, the
projected power is deliberately varied to check against model prediction, and
it is normally only conducted in the zenith without any requirement for
complex laser transfer optics for pointing mechanism, 2. to simplify the
calculation process. Only when discussing the annual mean on-sky performance
in different directions, we did the simulation for all five sites with the
full power of the laser at 20W for all elevation/azimuth angle combinations.

###  2.4 Temperature and atmospheric molecular density

Atmosphere molecular density and temperature in the mesosphere will influence
the collision possibility between sodium atoms and atmospheric molecules/atoms
and then affect the coupling efficiency and the photon return number
(Kibblewhite 2000, 2008). The values of the two parameters are related with
the location of a site and can be obtained by querying MSISE-90 model from
CCMC.

Site | Latitude | Longitude | Altitude |  Declination (+E $|$ -W) |  Inclination (+D $|$ -U) |  Strength of Local geomagnetic field (nT)
---|---|---|---|---|---|---
Xinglong | 40.39°N | 117.58°E | 960m | -7.05° | 59.33° | $5.20\times 10^{4}$
Lenghu | 38.45°N | 93.20°E | 2800m | 3.58° | 57.99° | $5.22\times 10^{4}$
Muztagh Ata | 38.28°N | 74.90°E | 4526m | 3.56° | 57.99° | $5.03\times 10^{4}$
Ali | 32.31°N | 80.05°E | 5100m | 1.61° | 50.33° | $4.80\times 10^{4}$
Gaomeigu | 26.71°N | 100.03°E | 3200m | -1.11° | 41.59° | $4.61\times 10^{4}$
Table 3: Parameters related to sites

##  3 RESULTS

###  3.1 Coupling efficiency variation under circular polarization at the
zenith

The coupling efficiency of circularly polarized light is higher than that of
linearly polarized light due to the larger scattering cross section of the
circular polarization. For this reason, all LGS lasers are circularly
polarized. In the simulation we take circular polarization as a default to
simulate the performance of the TIPC laser.

Figure 4: Coupling efficiency against projected power density on sky (92km) in
zenith direction

The comparison of the coupling efficiencies at the zenith with different power
densities is shown in Figure 4. The range of power density on sky (X-axis) is
from 0.01 to 1000 $W/m^{2}$ and the Y-axis is the coupling efficiency. The
curves of different colors represent the variations of coupling efficiency
against projected power densities at different sites. Coupling efficiencies at
Xinglong, Lenghu and Muztagh Ata are almost identical, which are higher than
Ali and Gaomeigu. This is mainly due to the influence of geomagnetic field.
When the direction of laser and geomagnetic field vector is parallel, lamor
precession caused by geomagnetic field effect is the least; when the two is
vertical, the influence of lamor precession is the strongest. From Table 3, we
can see that the angles of geomagnetic field and the zenith at the former
three sites are almost the same, and the geomagnetic field strength is
similar. Both of them are higher than that of the latter, which explains why
the coupling efficiency of Ali and Gaomeigu at zenith is lower.

###  3.2 Photon return efficiency in zenith

In the zenith direction, the sodium column density in the mesosphere will
change with month and latitude of the site, and the coupling efficiency will
change with the location of sites and power density. Therefore, we use the
product of the sodium column density and the coupling efficiency to
investigate the change of the photon return number with the power density in
different months at different sites. We temporarily call the product of the
two as the photon return efficiency in zenith.

Such photon return efficiencies from January to December are shown in Figure
5. Figure 5a shows the photon return efficiency variation with the power
density on sky in January. The curves with different colors represent the
efficiencies for different sites, which are the same as in Figure 4. Figure 5
b1-b12 shows the change of photon return efficiency from January to December
and the ranges of value for Y axis in these figures are the same for easy
comparison. We can summarize Figure 5 as follows:

  * •

In different months, the photon return efficiency of a site varies greatly.
For example, the maximum photon return efficiency of Xinglong in November is
129% higher than that in June: in June, the photon return efficiency is
$4.65\times 10^{15}$ (photons/s/W), but in November is $10.65\times 10^{15}$
(photons/s/W). This is mainly caused by the large change of sodium column
density with the month;

  * •

In spring and summer (from Mar to August), the photon return efficiencies of
different sites are similar and at a low level. While in autumn and winter
(from September to February), there is large differences between the photon
return efficiency of five sites. The higher the latitude of site, the higher
the photon return efficiency. There are two main reasons for this phenomenon:
one is that the coupling efficiency increases with latitude of site due to the
influence of geomagnetic field; the other is that the sodium column density
keeps at high level in autumn and winter.

Figure 5: Photon return efficiency variation against power density on sky
(92km) from Jan. to Dec (b1-b12).

###  3.3 Returned photon flux

The returned photon flux determines how much photons could reach the
telescope, which is exactly what astronomers are most concerned about. The
relationship between the coupling efficiency and the photon flux can be
described by equation 2 from Holzlöhner’s paper (Holzlöhner et al. 2010),
where Sce is the coupling efficiency, $\Phi$ is the photon flux received by
the telescope (unit photons/s/$m^{2}$), L is the vertical distance from the
primary of the telescope to the sodium layer, P is the laser power, T is the
atmospheric transparency, X = 1/($sec(\theta)$) and $\theta$ is the laser
launch zenith angle, $C_{Na}$ is the sodium column density. Assuming that the
atmospheric transparency at zenith is 90%, the seeing is 1.5” and beam quality
($M^{2}$) of TIPC laser is 1.4, the diameter of laser spot size at 92km on sky
is about 0.737m(1.652”), and the mean power density at zenith is about 42.2
W/$m^{2}$. The coupling efficiency can be obtained from the curve of the
coupling efficiency against power density(Figure 4). According to Equation 2,
the photon flux in zenith for different months can be calculated, as shown in
Figure 6 (left). In general, the photon flux is higher in winter and lower in
summer, and the flux variations at Xinglong, Lenghu and Muztagh Ata are
identical and higher than that of Ali and Gaomeigu, which is more pronounced
in winter. This is mainly due to the influence of the seasonal distribution of
the sodium layer column density. According to the relationship of photon flux
and magnitude of an object (Astronomical Magnitude Systems:
https://www.cfa.harvard.edu/ dfabricant/huchra/ay145/mags.html), we calculated
the corresponding Johnson V-band magnitude variation of the sodium guide star
with month, shown in Figure 6 (right). The highest V magnitude is up to 6.56
in November at Xinglong and the lowest is about 7.50 in June at Gaomeigu.

$\displaystyle\Phi=\frac{Sce\cdot P\cdot(T)^{2X}\cdot C_{Na}\cdot X}{L^{2}}$
(2)

Figure 6: Photon flux variation with month (left) and corresponding seasonal
distribution of V-band magnitude of sodium laser guide star (right)

###  3.4 Coupling efficiency variation with direction of the laser

As mentioned in the previous section, during an astronomical observation, the
direction of laser guide star would change with the pointing direction of a
telescope. Thus, the angle between the laser direction and the local
geomagnetic field vector will change, and this will cause continuous change of
coupling efficiency because of the varying influence of the geomagnetic field.
In addition, for slanted pointing directions, the increasing airmass will
reduce both the projected laser power and the return photon flux. Therefore,
the coupling efficiency has a peculiar spatial distribution on sky with
different pointing directions. In this case, we did our simulation with a
constant 20 W laser output power. For every laser pointing direction, we
calculated the angle between the laser and the geomagnetic field vector, and
applied this value into our simulation along with projected power density and
size of spot due to airmass and slanted projection. The atmospheric
transparency for every site is assumed to be 90%. The calculated on-sky
coupling efficiency results are shown in Figure 7.

In each picture of Figure 7, the X-axis represents the change of the azimuth
angle, with the true north at 0°and the east at 90°. The Y-axis represents the
change of the elevation angle, with the horizontal direction at 0°and the
zenith at 90°. The variation of the coupling efficiency with laser direction
at can be summarized as follows:

  * •

The coupling efficiency is almost symmetric along 180°of azimuth angle. When
the azimuth angle changes from north to south (0°-180°), the coupling
efficiency increases gradually and then decreases when the azimuth angle
changes from south to north (180°-360°);

  * •

The coupling efficiency increases first and then decreases with the increase
of elevation angle , but the elevation angles of the maximum coupling
efficiency are different among the sites;

  * •

In Xinglong, Lenghu and Muztagh Ata, the coupling efficiency reaches maximums
of 325.3, 327.2 and 325.5 (photons/s/W/(atoms/$m^{2}$)) at elevation angle of
60°and azimuth angle of about 180°; In Ali, the coupling efficiency reaches a
maximum of 326.1 (photons/s/W/(atoms/$m^{2}$)) at elevation angle of 50°and
azimuth angle of about 180°. In Gaomeigu, the coupling efficiency reaches a
maximum of 328.7 (photons/s/W/(atoms/$m^{2}$)) at elevation angle of 40°and
azimuth angle of about 180°. Statistics of the results are summarized in Table
4. The mean coupling efficiency is the average in all direction, and the
sodium column density is adopted the average of 12 months when we calculate
the corresponding V magnitude of mean coupling efficiency.

Site |  The elevation angle of the maximum coupling efficiency |  Mean coupling efficiency * (photons/s/W/ (atoms/$m^{2}$)) |  V magnitude of the mean coupling efficiency |  Maximum coupling efficiency (photons/s/W/ (atoms/$m^{2}$)) |  Minimum coupling efficiency (photons/s/W/ (atoms/$m^{2}$)) |  Elevation and azimuth angle of minimum coupling efficiency
---|---|---|---|---|---|---
Xinglong | 60° | 218.3 | 6.904 | 325.3 | 177.2 | (El:30°Az:330°)
Lenghu | 60° | 217.5 | 6.908 | 327.2 | 177.4 | (El:30°Az:40°)
Muztagh Ata | 60° | 218.9 | 6.901 | 325.5 | 178.7 | (El:30°Az:340°)
Ali | 50° | 216.1 | 6.915 | 326.1 | 181.4 | (El:30°Az:50°)
Gaomeigu | 40° | 212.1 | 6.935 | 328.7 | 182.5 | (El:30°Az:50°)
Table 4: Numerical information about coupling efficiency (* mean coupling
efficiency in all direction)

Xinglong

Lenghu

Muztagh Ata

Ali

Gaomeigu

Figure 7: Coupling efficiency distribution across the sky

##  4 Conclusion

In this paper, we have studied the performance of sodium laser guide stars
generated by a 20W class quasi-CW pulsed laser at five typical sites of China.
From Table 4, the overall annual mean performances of sodium guide star at
these sites are very similar, and the similarity can be summarized as due to
that the brightness of the guide star is strongly relying on the seasonal
variation of sodium column density and laser pointing direction rather than
geographical locations. However, for AO design concerns, coupling efficiency’s
variation in different months with different pointing angles are more
important because it determines the performance that the AO system could
function in specific night of year at a particular angle. The seasonal photon
flux distribution of LGS shows higher value in winter and lower in summer. The
photon return in summer is only about 44% (at Xinglong site) - 53% (at
Gaomeigu site) of that in winter (see Figure 6), which renders a difference of
0.882-0.685 in V magnitude. For different laser pointing direction, the
coupling efficiency changes most significantly with the azimuth angle at
elevation angle from 40°to 60°. In addition, the variation of coupling
efficiency with elevation angle is also prominent, which shows that the
elevation angle for the maximum coupling efficiency is higher with increasing
latitude as shown in Figure 7 and Table 4. With both effects combined, the
variation of sodium LGS’s brightness is not negligible in AO system’s design
process. We hope that this work could shed some light for future applications
of sodium laser guide star.

##  5 Acknowledgement

The research is supported by the Operation, Maintenance and Upgrading Fund for
Astronomical Telescopes and Facility Instruments, budgeted from the Ministry
of Finance of China (MOF) and administrated by the Chinese Academy of Sciences
(CAS).

## References

  * Bian et al. (2020) Bian, Q., Bo, Y., Zuo, J. W., et al., 2020, OptL, 45, 7, 1818
  * Bian et al. (2016a) Bian, Q., Bo, Y., Zuo, J. W., et al., 2016, OptL, 41, 8, 1732
  * Bian et al. (2016b) Bian, Q., Zuo, J. W., Guo, C., et al., 2016, LaPhy, 26, 095005
  * Bian et al. (2017) Bian, Q., Zong, Q. S., Chang, J. Q., et al., 2017, IPTL, 29, 2095
  * Calia et al. (2014) Calia, D. B. , Hackenberg, W. , Holzlöhner, R., et al., 2014, AdOT, 3, 3, 345
  * Chin et al. (2016) Chin, J. C. Y., Wizinowich, P., Wetherell, E., et al. 2016, Proc.SPIE, 9909, 99090S
  * Feng et al. (2016) Feng, L., Shen, Z. X., Xue, S. J., et al., 2016, RAA(Research in Astronomy and Astrophysics), 16, 144
  * Feng et al. (2015) Feng, L., Kibblewhite, E., Jin, K., et al., 2015, Proc. SPIE, 9678, 96781B
  * Feng et al. (2020) Feng, L., Hao, J. X., Cao, Z. H., et al., 2020, RAA(Research in Astronomy and Astrophysics), 20, 80
  * Feng & Hao (2020) Feng, L., Hao, J. X. 2020, RAA(Research in Astronomy and Astrophysics), 20, 79
  * Fussen et al. (2004) Fussen, D., Vanhellemont, F., Bingen, C., et al., 2004, GeoRL, 31, L24110
  * Fussen et al. (2010) Fussen, D., Vanhellemont, F., Tétard, C., et al., 2010, ACP, 10, 9225
  * Hardy (1998) Hardy, J. W., 1998, Adaptive Optics for Astronomical Telescopes (New York: Oxford University Press)
  * Holzlöhner et al. (2010) Holzlöhner, R., Rochester, S. M., Calia, D. B., et al., 2010, A&A, 510
  * Jin et al. (2015) Jin, K., Wei, K., Feng, L., et al., 2015, PASP, 127, 749
  * Jin et al. (2014) Jin, K., Wei, K., Xie, S., et al., 2014, Proc. SPIE, 9148, 91483L
  * Kibblewhite (2000) Kibblewhite, E. J., 2000, in Laser Guide Star Adaptive Optics for Astronomy, Ageorges, N. & Dainty, C., (Kluwer Academic Publishers), 51
  * Kibblewhite (2008) Kibblewhite, E. J., 2008, Proc. SPIE 7015, 70150M
  * Langowski et al. (2017) Langowski, M. P., Savigny, C., Burrows, J. P., et al., 2017, AMT, 10, 2989
  * Li et al. (2021) Li, H. Y., Feng, L., Wang, J. L., et al., 2021, PASP, 133, 1021
  * Liu et al. (2012) Liu, L. Y., Yao, Y. Q., Vernin, J., et al., 2012, in Proc. of SPIE, 8444, Ground-based and Airborne Telescope IV, 844464
  * Liu et al. (2015) Liu, L. Y., Yao, Y. Q., Vernin, J., et al., 2015, Journal of Physics: Conference Series, 595, 012019
  * Liu et al. (2016) Liu, Y., Song, T. F., Zhang, X. F., et al., 2016, in Proc. of IAU, 320, Solar and Stellar Flares and their Effects on Planets, 447
  * Moussaoui et al. (2008) Moussaoui, N., Holzlöhner, R., Hackenberg, W., et al., 2008, Proc. SPIE, 7015, 70152W
  * Moussaoui et al. (2009) Moussaoui, N., Holzlöhner, R., Hackenberg, W., et al., 2009, A&A, 501, 793
  * Otarola et al. (2016) Otarola, A., Hickson, P., Gagne, R., et al., 2016, JAI, 5, 1
  * Qian et al. (2015) Qian, X., Yao, Y. Q., Wang, H. S., et al., 2015, PKAS, 30, 695
  * Wang et al. (2015) Wang, H. S., Yao, Y. Q., Liu, L. Y., et al., 2015, Journal of Physics: Conference Series, 595, 012037
  * Wei et al. (2012) Wei, K., Bo, Y., Xue, X.H., et al. 2012, Proc.SPIE, 8447, 84471R
  * Wu et al. (2016) Wu, N., Liu, Y., & Zhang, H. M. 2016, AcASn, 57, 729
  * Xu et al. (2020) Xu, J., Esamdin, A., Hao, J.X., et al., 2020, RAA(Research in Astronomy and Astrophysics), 20, 86
  * Yao et al. (2013) Yao, Y. Q, Wang, Y. P., Liu, L. Y., et al., 2013, in NARIT Conf. Ser., 1, the 11th Asian-Pacific Regional IAU Meeting 2011, ed. Komonjinda, S. (IAU), 1, http://conference.narit.or.th/ncs/APRIM2011_proceedings/S7/S7P05.pdf


# Visualizing the (Causal) Effect of a Continuous Variable on a Time-To-Event
Outcome

Robin Denz and Nina Timmesfeld

Ruhr-University Bochum
Department of Medical Informatics, Biometry and Epidemiology

###### Abstract

Visualization is a key aspect of communicating the results of any study aiming
to estimate causal effects. In studies with time-to-event outcomes, the most
popular visualization approach is depicting survival curves stratified by the
variable of interest. This approach cannot be used when the variable of
interest is continuous. Simple workarounds, such as categorizing the
continuous covariate and plotting survival curves for each category, can
result in misleading depictions of the main effects. Instead, we propose a new
graphic, the _survival area plot_ , to directly depict the survival
probability over time and as a function of a continuous covariate
simultaneously. This plot utilizes g-computation based on a suitable time-to-
event model to obtain the relevant estimates. Through the use of
g-computation, those estimates can be adjusted for confounding without
additional effort, allowing a causal interpretation under the standard causal
identifiability assumptions. If those assumptions are not met, the proposed
plot may still be used to depict non-causal associations. We illustrate and
compare the proposed graphics to simpler alternatives using data from a large
German observational study investigating the effect of the Ankle Brachial
Index on survival. To facilitate the usage of these plots, we additionally
developed the contsurvplot R-package which includes all methods discussed in
this paper.

_Keywords: continuous variables, time-to-event outcome, confounder-adjustment,
counterfactual survival probability_

##  1 Introduction

Visualization is a crucial aspect when reporting the main results of any
study. When the aim of the study is to estimate a causal effect of one
variable on a specific endpoint, a clear graphical display of the effect can
make the results more accessible for the reader [56]. With a time-to-event
endpoint, such as time until recurrence of a disease or death, the most
popular approach is to draw Kaplan-Meier survival curves, or confounder-
adjusted alternatives, for each group of interest [11]. This, however, only
works if there is a discrete number of groups. In reality, the variable of
interest may be continuous, making the use of Kaplan-Meier survival curves
impossible. A few examples for continuous covariates that are frequently of
interest are income, duration of exercise or tissue size [44, 51].

As a workaround, researchers often resort to either (1) performing artificial
categorization in both the main analysis and the graphical displays, or (2)
performing categorization only for the graphics while using the full
continuous variable in the main analysis. Dividing the variable into
artificial groups at arbitrary or “optimal” values [31, 18] has the advantage
that survival curves can be created. It is, however, problematic due to a loss
of information, the seemingly arbitrary choice of cut-off points and the
improbable assumption of homogeneity of risk inside the categories [36, 4, 3].
In case (1) the statistical power of the analysis might be reduced and or bias
might be introduced. In case (2), the survival curves might show a misleading
picture of the treatment effect, which may also not be consistent with the
results of the main analysis because different scales were used [31].

Graphics designed for displaying the effect of a continuous covariate on a
time-to-event outcome exist, but are rarely used in practice [26, 34, 14, 47,
48, 24]. Most of these approaches fail to represent the treatment effect over
time and as a function of the continuous variable simultaneously. Some rely on
summary statistics, such as landmark survival probabilities or restricted mean
survival times [47, 14], while others rely on rather complex estimation
strategies that may not use easily interpretable absolute measures of effect
[26, 34, 48, 24]. Furthermore, most of the proposed graphics cannot be
adjusted for confounders and can therefore only be used to display causal
effects in the absence of confounding. This is rarely the case when the
treatment is continuous, because continuous variables are usually not randomly
assigned.

In this article, we propose a new type of graphic, the _survival area plot_ ,
that does not suffer from these problems. The survival area plot directly
extends the standard Kaplan-Meier plot to the continuous case, by depicting
the survival probability as a function of time and a continuous covariate
through the use of a colour scaled area. We use g-computation [43] to obtain
the required probability estimates, allowing the plot to be adjusted for
confounders. Under the standard causal identifiability assumptions [37, 21],
the resulting survival probabilities can be interpreted as counterfactual
quantities, resulting in a display of the causal effect of the continuous
variable on the time-to-event outcome. To facilitate the use of this method in
practice, we also developed the R-package contsurvplot [10], containing an
implementation of the proposed and alternative methods in the R programming
language.

First, we give a formal description of the target estimand and the estimation
methodology. Next, we introduce the survival area plot and describe
alternative ways to visualize the g-computation based estimates. Afterwards,
examples of the graphics are presented using data from the German
Epidemiological Trial on Ankle Brachial Index [17]. Using these examples, we
compare the different approaches and justify our practical recommendations.

##  2 Proposed Method

###  2.1 Target Estimand

Before going further, it is crucial to define what “the causal effect of a
continuous variable on a time-to-event outcome” is. We will first discuss the
target estimand in the simpler case of a dichotomous variable. Let $Z$ denote
the variable of interest with $Z\in\\{0,1\\}$. Let $T$ be the time until the
occurrence of an event. In this setup, each individual has two potential
failure times: $T^{Z=0}$ being the time until the event occurs under
treatment-regime $Z=0$ and $T^{Z=1}$ being the failure time under $Z=1$. Our
goal is to estimate the population-level counterfactual survival function
under both $Z=0$ and $Z=1$, defined as:

$S_{z}(t)=E\left(I\left(T^{Z=z}>t\right)\right).$ (1)

This is the survival probability that would have been observed at time $t$ if
the value of $Z$ of all individuals in the population of interest had been set
to $z$ (_counterfactual survival probability_) [11]. If $Z$ was randomized and
the study was performed adequately, this is what a Kaplan-Meier estimator [25]
estimates. The difference between treatment-specific counterfactual quantities
($S_{1}(t)-S_{0}(t)$) or their ratio ($S_{1}(t)/S_{0}(t)$) may be used to
define average causal effects [50].

It is straightforward to generalize this notation to the continuous case [15].
Let $Z$ be a continuous variable that can take any real values
$Z\subseteq\mathbb{R}$. Instead of two potential failure times there is now an
uncountably infinite set of potential failure times $T^{Z=z}$ for each
individual, but equation 1 and its interpretation remain unchanged. Average
causal effects may still be defined using differences or ratios, as explained
in the e-appendix. Note that this counterfactual quantity is only meaningful
if $Z$ corresponds to a well-defined treatment, which may not always be the
case with continuous $Z$ [20, 13, 49].

###  2.2 Estimation using G-Computation

Suppose that setting $Z$ to $z$ corresponds to a well-defined intervention. If
$Z$ was randomly assigned, estimating $S_{z}(t)$ is “only” an issue of
interpolation. As mentioned above if $Z$ was dichotomous and appropriately
randomized, we could obtain an unbiased estimate of $S_{z}(t)$ using a
standard stratified Kaplan-Meier estimator [25]. This, however, does not work
with a continuous $Z$, because of its infinite amount of possible values. A
suitable model is necessary to estimate this function. A popular alternative
is the Cox proportional hazards regression model [9], given by:

$h(t)=h_{0}(t)\exp(\beta Z),$ (2)

where $\beta$ is the coefficient of the continuous variable and $h_{0}(t)$ is
the baseline-hazard function, which can be estimated from the data [5]. This
model can be used to predict the survival probability of an individual with a
specific value of $Z$ at a given point in time $t$ using:

$S(t|Z=z)=\exp\left(-\int_{0}^{t}h_{0}(t)e^{\beta Z}dt\right).$ (3)

If $Z$ was randomized and the required assumptions for the Cox model hold,
using a Cox model that includes only $Z$ as an independent variable and an
estimate of $h_{0}(t)$ is enough to obtain an unbiased estimate for
$S_{z}(t)$. Without randomization, we also need to take confounding into
account. This can be done by including a sufficient set of confounders in the
Cox model as additional independent variables (with an appropriate functional
form) and using this model to perform _regression standardization_ , also
known as _g-computation_ [43, 27]. Here, the value for $Z$ is set to $z$ for
each person in the sample. The survival probability at time $t$ is then
estimated for each individual person by using the model and the observed
values of the included confounders. The mean of these person-specific
estimates in a sample of size $n$ is an unbiased estimator of $S_{z}(t)$.
Formally, this can be defined as:

$\hat{S}_{z}(t)=\frac{1}{n}\sum_{i=1}^{n}S(t|Z=z,X=x_{i}),$ (4)

where $X$ is a vector of confounders and $x_{i}$ is the realization of this
vector for the $i$-th person. $S(t|Z=z,X=x_{i})$ corresponds to equation (3),
but including a set of confounders $X$. This method has been used extensively
when analyzing time-to-event data [30, 7]. Note that it is not necessary to
use a Cox model. Any model capable of making predictions of $S(t|Z=z,X=x_{i})$
may be used instead. Possible alternatives include additive hazards models
[1], parametric survival time models [23] or random survival forests [22].

Four fundamental assumptions have to be met to allow a causal interpretation
of the obtained values. Firstly, the _no interference assumption_ states that
the potential survival time of an individual is independent of the value of
$Z$ of other individuals [37]. Secondly, _counterfactual consistency_ posits
that the potential survival time would remain the same under $Z$ regardless of
whether $Z$ was set experimentally or not ($P(T=t|Z=z)=P(T^{Z=z}=t|Z=z)$) [20,
40]. Thirdly, the _conditional exchangeability assumption_ asserts that the
treatment groups are exchangeable, given a sufficient set of relevant
confounders, i.e.,
$T^{Z=z}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}Z|X$ [46].
Lastly, the _positivity assumption_ states that every individual has a non-
zero probability of receiving any value of $Z$ ($P(Z=z|X)>0$ for all
$Z\subseteq\mathbb{R}$). Although it is impossible to observe individuals for
every value of $Z$ due to the uncountably infinite number of possible values,
this _practical violation_ is generally not an issue because of the
interpolation performed by the time-to-event model [52, 27, 21].

In addition any model-specific assumptions, such as the proportional hazards
assumption when using a Cox model [9], also have to be met. This includes the
correct specification of the functional form of all included variables and the
proper treatment of censored observations. Throughout this article, we assume
random right-censoring. A discussion about the implications of dependent
censoring is given in the e-appendix. If any of these assumptions are
violated, the resulting estimates cannot be interpreted as counterfactual
survival probabilities. However, the graphics proposed in this article may
still be used to visualize the non-causal marginal association of the
continuous covariate and the time-to-event endpoint. This may be useful to
communicate results of prediction models, especially when using machine
learning methods that are otherwise hard to interpret [22, 42].

###  2.3 Computational Details

Our implementation [10] of g-computation works by first setting $Z$ to $z$ for
all individuals included in the sample. Afterwards, the time-to-event model is
used to obtain predictions of $S(t|Z=z,X=x_{i})$ for each individual at $t$.
Internally, the fast predictRisk function from the riskRegression R-package
[38] is used at this step. Finally, those predictions are averaged for each
considered $t$. If multiple values of $Z$ are supplied, this procedure is
simply repeated. This algorithm can be slow when considering many values of
$t$ and $Z$ simultaneously. Nevertheless, it still executes in reasonable time
($<1$min) for medium sized datasets ($n<5000$) on a regular computer. The
runtime of this procedure may be decreased by executing the computations for
different values of $Z$ on multiple cores simultaneously.

###  2.4 Visualization of the Results

Performing the procedure described above for multiple $t$ and several values
for $Z$ results in a three-dimensional survival surface, which can be used to
create multiple plots based on summary statistics. Most of these plot have
already been described in the literature [54, 47, 14, 24], but the authors
introducing them did not use g-computation to obtain the underlying estimates.
By using g-computation instead, the values represented by the graphs can be
interpreted as counterfactual quantities, if the previously mentioned
assumptions hold.

The simplest of these is a _landmark survival probability plot_ [54, 47]. To
create this plot, one or multiple $t$ are chosen. The survival probability at
$t$ is then visualized as a function of the continuous covariate. A _survival
time quantile plot_ follows the same structure, but replaces the landmark
survival probability with a chosen survival time quantile. The same type of
plot can also be created using the _restricted mean survival time_ [14].
Another alternative would be to plot standard survival curves, for specific
values of $Z$ (_value specific survival curves_) [14].

All of these techniques reduce the surface to two dimensions. Although this
can be an effective strategy in some situations, we argue that the most
appropriate visualization technique is to represent the whole survival surface
instead. It could be argued that the most natural way to present a three-
dimensional surface would be a three-dimensional plot [48], but this is not
what we advocate for. Many researchers have pointed out the shortcomings of
non-interactive three-dimensional plots before [45, 53]. Those include the
difficulty in reading off specific values of the surface and the seemingly
inconsistent look of the graphic when viewed from different angles [53].

Instead, the surface can be represented in a two-dimensional plot using color
scales. We advocate for two options: _survival contour plots_ first introduced
by [24] and our own proposal, _survival area plots_. Survival contour plots
display the time on the $x$-axis and the continuous covariate on the $y$-axis.
The survival probability is categorized and the corresponding area is colored
accordingly. The original version proposed by [24] used a kernel-density based
estimation strategy, which, contrary to our g-computation based approach, does
not allow confounder-adjustment.

Survival area plots keep the structure of Kaplan-Meier plots with $t$ on the
$x$-axis and the survival probability on the $y$-axis. Instead of single
survival curves for some categories, an area is drawn with the filled color
corresponding to values of $Z$. This can be a continuously colored area or
divided into bins similar to the survival contour plot. Essentially, survival
area plots depict all possible _value specific survival curves_ in a given
range in the same graphic. Since those change continuously over changing
values of $Z$, they create an area instead of discernible lines. Therefore,
every infinitely thin part of this area is equal to one specific estimate of
$S_{z}(t)$ over the depicted time range. Because those individual parts can be
interpreted as the survival probability that would have been observed if $Z$
was set to $z$ in the target population, the area depicts how this probability
changes with changing values of $Z$.

##  3 Illustrative Examples

To illustrate our proposed methodology, we use data from the German
Epidemiological Trial on Ankle Brachial Index (getABI) [17]. The getABI study
is a prospective observational study including 6880 primary care patients aged
65 or older who were followed for up to 8 years. An initial examination was
performed for every patient in the beginning of the study in 2001. In this
examination, both the Ankle-Brachial-Index (ABI) and other clinically relevant
factors, such as past cardiovascular events and socio-economic factors, were
collected. The ABI was of special interest in this study, because it is a
frequently used marker to diagnose peripheral artery disease (PAD). It is
calculated as “the ratio of the higher of the two systolic pressures (tibial
posterior and anterior artery) above the ankle to the average of the right and
left brachial artery pressures” [12].

We want to estimate the fraction of patients that would have survived between
0 and 8 years if the ABI of every member in the target population had been set
to values between 0.35 and 1.5. Because the ABI cannot be “set” to a specific
value in reality, this is a hypothetical intervention. It is also somewhat
imprecise, because it does not refer to any duration of time in which the ABI
stays at that level. Although some imprecision in the definition is inevitable
[20, 49] it is still important consider these limitations when interpreting
the results. Our target population consists of all German primary care
patients aged 65 or older with a prior cardiovascular event at baseline. In
accord to previous studies, we also only consider patients with an ABI of
$<=1.5$ at baseline, because an ABI $>1.5$ indicates incompressible arteries
[33, 32]. Applying these restrictions to the getABI data results in a sample
size of 508 patients with 143 events.

In practice, the ABI is often categorized, indicating different levels of PAD.
An ABI $<0.9$ is usually considered as “low” and an indication of asymptomatic
PAD [12, 35, 29]. Some authors consider an ABI between $0.9$ and $1.3$ as
“normal” [35], while others extend this category up to an ABI $1.4$ or $1.5$
[12, 29]. These categorizations are frequently used when analyzing the effect
of the ABI on the survival time, which, as stated previously, can lead to a
loss of statistical power and might even result in misleading conclusions
[31].

To illustrate the difference between the classical and our proposed methods,
we first fitted five Cox models. All of these models use death of all causes
as response variable, but utilize different specifications of the ABI. The
first model includes the ABI as a categorical variable and the relevant
confounders age, sex and smoking status. The second one contains only the ABI
as independent continuous covariate. The third one adds the relevant
confounders to model two. Model (4) and (5) are identical to models (2) and
(3), with the only difference being that ABI was modeled using B-Splines [41],
allowing a non-linear fit.

Table 1: Hazard Ratios and associated 95% confidence intervals of different
Cox proportional hazards regression models describing the relationship of the
ABI on the survival time.

| Dependent variable:
---|---
| Death of all causes
| (1) | (2) | (3) | (4) | (5)
Age | 1.091 |  | 1.091 |  | 1.089
| (1.056, 1.127) |  | (1.056, 1.127) |  | (1.054, 1.125)
Sex | 0.988 |  | 1.009 |  | 0.983
| (0.661, 1.477) |  | (0.676, 1.507) |  | (0.657, 1.469)
Ever Smoker | 2.202 |  | 2.208 |  | 2.225
| (1.410, 3.440) |  | (1.411, 3.454) |  | (1.424, 3.476)
ABI (0.5,0.7] | 0.705 |  |  |  |
| (0.298, 1.667) |  |  |  |
ABI (0.7,0.9] | 0.524 |  |  |  |
| (0.243, 1.131) |  |  |  |
ABI (0.9,1.1] | 0.366 |  |  |  |
| (0.173, 0.773) |  |  |  |
ABI (1.1,1.5] | 0.364 |  |  |  |
| (0.167, 0.793) |  |  |  |
ABI (per 0.1 points) |  | 0.867 | 0.894 |  |
|  | (0.802, 0.937) | (0.826, 0.967) |  |
bs(ABI, df = 3)1 |  |  |  | 0.402 | 0.286
|  |  |  | (0.037, 4.343) | (0.026, 3.148)
bs(ABI, df = 3)2 |  |  |  | 0.098 | 0.152
|  |  |  | (0.022, 0.430) | (0.036, 0.653)
bs(ABI, df = 3)3 |  |  |  | 0.380 | 0.378
|  |  |  | (0.066, 2.176) | (0.069, 2.082)
Observations | 508 | 508 | 508 | 508 | 508
R2 | 0.101 | 0.024 | 0.095 | 0.032 | 0.102
Max. Possible R2 | 0.966 | 0.966 | 0.966 | 0.966 | 0.966
Log Likelihood |  $-$831.596 |  $-$852.444 |  $-$833.311 |  $-$850.397 |  $-$831.343

The hazard ratios of all models are displayed in table 1. Note that these
hazard ratios can generally not be endowed with a causal interpretation [19,
2]. However, it is still apparent that the ABI is associated with the survival
time, regardless of the model specification used. One can also see that once
age, sex, and smoking status are controlled for, the association becomes
smaller. While we can infer these trends from the hazard ratios, it is still
difficult to accurately interpret the results, especially when B-Splines were
used to model the ABI. Graphical displays of the survival probabilities are
necessary to enhance understanding here.

Figure 1 displays survival curves stratified by the ABI categories defined by
[12]. The left panel shows simple Kaplan-Meier estimates, stratified by the
ABI categories, while the right panel shows survival curves adjusted for sex,
age and smoking status. To adjust these curves, we used model (1) to perform
g-computation. A clear dose-response effect can be seen, regardless of whether
the curves were adjusted for confounders or not. The category with the lowest
ABI values has the lowest survival probability, while the categories with
higher ABI values have a higher survival probability. In the adjusted graph,
there is no visible difference between the survival probability of the two
upper categories. Looking only at Figure 1 or model 1, it remains unclear
whether the survival probability remains constant with increasing ABI values
after a certain point or if there is a non-linear relationship. It also
implies that the risk associated with an ABI of 0.7 is equal to the risk
associated with an ABI of 0.9, due to the assumed homogeneity inside each
category [4].

Figure 1: Survival curves stratified by categories of the ABI. Results in the
left panel were obtained using a stratified Kaplan-Meier estimator while the
results in the right panel were obtained using g-computation and model (1)
displayed in table 1. The survival curves for the highest two categories are
displayed as transparent in the right panel, to make it more obvious that
those are two curves.

Instead of using the categorical version of the ABI variable, we can also use
model (2) and (3) to create survival area plots. Figure 2 displays two kinds
of these plots. On the left side a continuous color scale is used to represent
the ABI values, while the right side shows the same plots using a discrete
color scale. The version on the left highlights the continuous nature of the
causal effect. It is, however, difficult to discern survival probabilities for
specific values of the ABI. The right version explicitly shows the range in
which the survival probability for the respective categories is in. We argue
that both versions give a more accurate depiction of the underlying treatment
effect than the simple categorized survival curves, because they correctly
depict the change of the survival probability with changing ABI values instead
of showing homogeneous probabilities in some categories.

Figure 2: Continuously and discretely colored survival area plots displaying
different estimates of the causal effect of the ABI on survival, created using
different Cox proportional hazards regression models and g-computation.
A: Using model 1 (unadjusted), continuous color shading; B: Using model 1
(unadjusted), discrete color shading; C: Using model 2 (adjusted), continuous
color shading; D: Using model 2 (adjusted), discrete color shading.

In models (2) and (3), we assumed that the ABI has a linear effect on the log
scale. The categorized plot shows some evidence that this assumption may be
wrong when the ABI values are high ($>1.1$), which is why we fitted models (4)
and (5) where we used B-Splines to model the ABI. Unfortunately, we cannot use
these models to create standard survival area plots, because they cannot
depict non-monotonic effects. If there is a curved relationship, the
corresponding areas would lie on top of each other, rendering the graphic
unusable. The only way to still use survival area plots in this case is to
divide the plot into multiple facets, where each facet contains only a part of
the entire surface. We discuss this in detail in the e-appendix.

Instead, survival contour plots can be used. Figure 3 shows the corresponding
plots for models (2) to (5). As before, there are only small differences
between the adjusted and the unadjusted analysis. However, a clear difference
can be seen when comparing the survival contour plots created using the
linearity assumption with the ones that were created using B-spline models.
There seems to be a non-linear relationship, where the survival probability
decreases again with ABI values higher than approximately 1. It is unlikely
that this is an artifact due to the population size, because there are 290
patients corresponding to 68 events, with an ABI $>1$.

Figure 3: Survival contour plots displaying the causal effect of the ABI on
survival, created using four different Cox proportional hazards regression
models and g-computation.

##  4 Discussion

In this article, we proposed the _survival area plot_ , which is a novel
method to visualize the effect of a continuous variable on a time-to-event
outcome. Additionally, we extended existing methods, such as _survival contour
plots_ [24] to allow confounder-adjustment through the use of g-computation.
We illustrated these methods using data from the German epidemiological trial
on Ankle Brachical Index [17] and compared them with the standard method of
categorized survival curves. Both survival area plots and survival contour
plots correctly display the estimated counterfactual survival probability as a
function of the continuous variable and time simultaneously. They avoid the
problems associated with categorization of continuous variables, such as loss
of statistical power, bias and the assumption of homogeneity in intervals [3,
4]. The survival area plot has the advantage that it keeps the well known
Kaplan-Meier plot structure and might therefore be easier to understand. It
is, however, difficult to draw this kind of plot if the effect of the
continuous variable is not monotonic. Survival contour plots do not share this
disadvantage, allowing any kind of non-linear effect.

Using g-computation to obtain the required estimates is both a strength and a
weakness for multiple reasons. First, it can be considered a strength because
it is widely used, relatively simple to understand and implement. For most
clinical researcher it is also straightforward to include into their data
analysis workflow, as they are already using regression models such as the Cox
model. Additionally, it has been shown to be one of the most efficient methods
when used properly [8, 11]. However, the method requires that a valid model
can be fit to the data, which might be difficult or even impossible in some
situations.

Like other methods used for causal inference, untestable assumptions have to
be made in order to allow a counterfactual interpretation of the estimates.
Among those is the assumption that the variable of interest constitutes a
“well-defined treatment”, which may be violated for some continuous variables
such as the Body-Mass-Index. [49, 13, 16] This is, however, not always the
case. Many variables, such as monthly income or duration of exercise, are both
directly manipulable and clearly well-defined treatments. A critical case-by-
case assessment of the relevant assumptions is necessary to allow a causal
interpretation of the resulting estimates and subsequently of our proposed
graphics. If any of the assumptions seem unrealistic, the obtained estimates
can only be considered to be marginal associations.

Another limitation of area-based plots is that there is no straightforward way
to display uncertainty in them. In standard Kaplan-Meier plots, confidence
intervals can be added to give an estimate of the statistical uncertainty of
the estimates. This is impossible in graphics based an colored areas. Plots
based on summary statistics might be preferable when the uncertainty of the
estimates is high. Nevertheless, we still recommend basing these plots on
g-computation when possible. This can be done using our contsurvplot R-package
[10], which contains implementations for all plots mentioned in this article.

We only considered a time-to-event scenario with a single type of event in
this article. Left-truncation or interval-censoring can be important
considerations in practice [28]. It is straightforward to use our proposed
methods even in these or other scenarios. All that has to be done is to fit a
suitable model which can be used to perform the required g-computation step.
That model may include any of the features mentioned. Extending the proposed
methodology to the context of time-varying treatments or settings with
competing events of interest is more complex, because it would require a
different target estimand [55]. Further work is necessary to extend the
proposed methods to those scenarios.

Furthermore, g-computation is only one of many methods to estimate
counterfactual survival probabilities [11]. It is a single robust method,
which means that it relies on the outcome model to obtain unbiased results. To
create the plots described in this paper, it is not necessary to rely on this
particular estimation method. Any method may be used, as long as it allows
estimation of the defined target estimand. Different authors recommend using
doubly-robust methods, such as Augmented-Inverse-Probability of Treatment
Weighting [39] or Targeted Maximum Likelihood Estimation [6] for similarly
defined estimands. Those methods rely on both an outcome model and a
treatment-assignment model and produce unbiased results whenever at least one
of these models is correctly specified. Using these methods would be a
promising alternative to g-computation. However, this is currently not
possible, because they do not support continuous exposures. Further research
is necessary to extent those estimators to this special case.

## Acknowledgments

We would like to thank Prof. Dr. Andreas Stang, Marianne Tokic and Prof. Dr.
Hans J. Trampisch for their helpful comments and suggestions. Additionally, we
want to thank the anonymous reviewers for their extensive comments and
suggestions which greatly improved the manuscript and the associated
R-package.

## References

  * [1] Odd Aalen “A Model for Nonparametric Regression Analysis of Counting Processes” In _Mathematical Statistics and Probability Theory: Lecture Notes in Statistics_ New York: Springer, 1980
  * [2] Odd O. Aalen, Richard J. Cook and Kjetil Røysland “Does Cox Analysis of a Randomized Survival Study Yield a Causal Treatment Effect?” In _Lifetime Data Analysis_ 21, 2015, pp. 579–593
  * [3] Douglas G. Altman “Categorizing Continuous Variables” In _Wiley StatsRef: Statistics Reference Online_ John Wiley & Sons, 2014
  * [4] Caroline Bennette and Andrew Vickers “Against Quantiles: Categorization of Continuous Variables in Epidemiologic Research, and Its Discontents” In _BMC Medical Research Methodology_ 12.21, 2012
  * [5] N. Breslow “Discussion of the Paper by D. R. Cox” In _Journal of the Royal Statistical Society: Series B_ 34.2, 1972, pp. 216–217
  * [6] Weixin Cai and Mark J. van der Laan “One-Step Targeted Maximum Likelihood Estimation for Time-To-Event Outcomes” In _Biometrics_ 76, 2020, pp. 722–733
  * [7] I-Ming Chang, Rebecca Gelman and Marcello Pagano “Corrected Group Prognostic Curves and Summary Statistics” In _Journal of Chronic Diseases_ 35, 1982, pp. 669–674
  * [8] Arthur Chatton, Florent Le Borgne, Clémence Leyrat, Florence Gillaizeau, Chloé Rousseau, Laetitia Barbin, David Laplaud, Maxime Léger, Bruno Giraudeau and Yohann Foucher “G-Computation, Propensity Score-Based Methods, and Targeted Maximum Likelihood Estimator for Causal Inference with Different Covariate Sets: A Comparative Simulation Study” In _Scientific Reports_ 10.9219, 2020
  * [9] D. R. Cox “Regression Models and Life-Tables” In _Journal of the Royal Statistical Society: Series B (Methodological)_ 34.2, 1972, pp. 187–220
  * [10] Robin Denz “contsurvplot: Visualize the Effect of a Continuous Varibale on a Time-To-Event Outcome” https://cran.r-project.org/package=contsurvplot, R package version 0.1.0, 2022
  * [11] Robin Denz, Renate Klaaßen-Mielke and Nina Timmesfeld “A Comparison of Different Methods to Adjust Survival Curves for Confounders” In _Statistics in Medicine_ InPrint, 2023
  * [12] Curt Diehm, Alexander Schuster, Jens R. Allenberg, Harald Darius, Roman Haberl, Stefan Lange, David Pittrow, Berndt Stritzky, Gerhart Tepohl and Hans-Joachim Trampisch “High Prevalence of Peripheral Arterial Disease and Co-Morbidity in 6880 Primary Care Patients: Cross-Sectional Study” In _Atherosclerosis_ 172.1, 2004, pp. 95–105
  * [13] Zili Dong “Well-Defined Interventions and Causal Variable Choice” In _Philosophy of Science_ InPrint, 2022
  * [14] Keving H. Eng, Emily Schiller and Kayle Morrell “On Representing the Prognostic Value of Continuous Gene Expression Biomarkers with the Restricted Mean Survival Curve” In _Oncotarget_ 6.34, 2015, pp. 36308–36318
  * [15] Douglas Galagate “Causal Inference with a Continuous Treatment and Outcome: Alternative Estimators for Parametric Dose-Response Functions with Applications”, 2
  * [16] Mekdes K. Gebremariam, Roch A. Nianogo, Nanna Lien, Mona Bjelland, Knut-Inge Klepp, Ingunn H. Bergh, Yngvar Ommundsen and Onyebuchi A. Arah “The Impact of Hypothetical Interventions on Adiposity in Adolescence” In _Scientific Reports_ 11.11216, 2021
  * [17] getABI Study Group “getABI: German Epidemiological Trial on Ankle Brachial Index for Elderly Patients in Family Practice to Dedect Peripheral Arterial Disease, Significant Marker for High Mortality” In _VASA_ 31.4, 2002, pp. 241–248
  * [18] Alberto Giannoni, Resham Baruah, Tora Leong, Michaela B. Rehman, Luigi Emilio Pastormerlo, Frank E. Harrell, Andrew J. S. Coats and Darrel P. Francis “Do Optimal Prognostic Thresholds in Continuous Physiological Variables Really Exist? Analysis of Origin of Apparent Thresholds, with Systematic Review for Peak Oxygen Consumption, Ejection Fraction and BNP” In _PLoS One_ 9.1, 2014
  * [19] Miguel A. Hernán “The Hazards of Hazard Ratios” In _Epidemiology_ 21.1, 2010, pp. 13–15
  * [20] Miguel A. Hernán “Does Water Kill? A Call for Less Casual Causal Inferences” In _Annual Journal of Epidemiology_ 26.10, 2016, pp. 674–680
  * [21] Miguel A. Hernán and James M. Robins “Causal Inference: What If” CRC Press, 2020
  * [22] Hemant Ishwaran, Udaya B. Kogalur, Eugene H. Blackstone and Michael S. Lauer “Random Survival Forests” In _The Annals of Applied Statistics_ 2.3, 2008, pp. 841–860
  * [23] Christopher H. Jackson “flexsurv: A Platform for Parametric Survival Modeling in R” In _Journal of Statistical Software_ 70.8, 2016
  * [24] Richard J. Jackson and Trevor F. Cox “Kernel Hazard Estimation for Visualisation of the Effect of a Continuous Covariate on Time-To-Event Endpoints” In _Pharamaceutical Statistics_ 21.3, 2022, pp. 514–524
  * [25] E. L. Kaplan and Paul Meier “Nonparametric Estimation from Incomplete Observations” In _Journal of the American Statistical Association_ 53.282, 1958, pp. 457–481
  * [26] Juha Karvanen and Frank E. Harrell Jr. “Visualizing Covariates in Proportional Hazards Model” In _Statistics in Medicine_ 28, 2009, pp. 1957–1966
  * [27] Alexander P. Keil, Jessie K. Edwards, David R. Richardson, Ashley I. Naimi and Stephen R. Cole “The Parametric G-Formula for Time-To-Event Data: Towards Intuition with a Worked Example” In _Epidemiology_ 25.6, 2014, pp. 889–897
  * [28] David G. Kleinbaum and Mitchel Klein “Survival Analysis: A Self-Learning Text” New York: Springer, 2012
  * [29] Jarosław Królczyk, Anna Skalska, Karolina Piotrowic, Małgorzata Mossakowska, Tomasz Grodzicki and Jerzy Gąsowksi “Disparate Effects of Ankle-Brachial Index on Mortality in the ’very old’ and ’younger old’ populations: The PolSenior Survey” In _Heart and Vessels_ 37, 2022, pp. 665–672
  * [30] Robert W. Makuch “Adjusted Survival Curve Estimation Using Covariates” In _Journal of Chronic Diseases_ 35.6, 1982, pp. 437–443
  * [31] Madhu Mazumdar and Jill R. Glassman “Categorizing a Prognostic Variable: Review of Methods, Code for Easy Implementation and Applications to Decision-Making about Cancer Treatments” In _Statistics in Medicine_ 19, 2000, pp. 113–132
  * [32] Mary McGrae McDermott, Philip Greenland, Kiang Liu, Jack M. Guralnik, Lilian Celic, Michael H. Criqui, Cheeling Chang, Gary J. Martin, Joseph Schneider, William H. Pearce, Lloyd M. Taylor and Elizabeth Clark “The Ankle Brachial Index is Associated with Leg Function and Physical Activity: The Walking and Leg Circulation Study” In _Annals of Internal Medicine_ 136.12, 2002, pp. 873–883
  * [33] Wouter T. Meijer, Arno W. Hoes, Dominique Rutgers, Michiel L. Bots, Albert Hofman and Diederick E. Grobbee “Peripheral Arterial Disease in the Elderly: The Rotterdam Study” In _Arteriosclerosis, Thrombosis, and Vascular Biology_ 18.2, 1998, pp. 185–192
  * [34] Luís Meira-Machado, Carmen Cardarso-Suárez, Francisco Gude and Artur Araújo “smoothHR: An R Package for Pointwise Nonparametric Estimation of Hazard Ratio Curves of Continuous Predictors” In _Computational and Mathematical Methods in Medicine_ , 2013
  * [35] Jair Baptista Miguel, Jorge Paulo Strogoff Matos and Jocemir Ronaldo Lugon “Ankle-Brachial Index as a Predictor of Mortality in Hemodialysis: A 5-Year Cohort Study” In _Arquivos Brasileiros de Cardiologia_ 108.3, 2017, pp. 204–211
  * [36] O. Naggara, J. Raymond, F. Guilbert, D. Roy, A. Weill and D. G. Altman “Analysis by Categorizing or Dichotomizing Continuous Variables is Inadvisable: An Example from the Natural History of Unruptered Aneurysms” In _American Journal of Neuroradiology_ 32.3, 2011, pp. 437–440
  * [37] Ashley I. Naimi and Jay S. Kaufman “Counterfactual Theory in Social Epidemiology: Reconciling Analysis and Action for the Social Determinants of Health” In _Current Epidemiology Reports_ 2, 2015, pp. 52–60
  * [38] Brice Ozenne, Anne Lyngholm Sørensen, Thomas Scheike, Christian Torp-Pedersen and Thomas Alexander Gerds “riskRegression: Predicting the Risk of an Event using Cox Regression Models” In _The R Journal_ 9.2, 2017, pp. 440–460
  * [39] Brice Maxime Hugues Ozenne, Thomas Harder Scheike and Laila Stærk “On the Estimation of Average Treatment Effects with Right-Censored Time to Event Outcome and Competing Risks” In _Biometrical Journal_ 62, 2020, pp. 751–763
  * [40] Judea Pearl “Does Obesity Shorten Life? Or is it the Soda? On Non-Manipulable Causes” In _Journal of Causal Inference_ 6.2, 2018
  * [41] Aris Perperoglou, Willi Sauerbrei, Michal Abrahamowic and Matthias Schmid “A Review of Spline Function Procedures in R” In _BMC Medical Research Methodology_ 19.46, 2019
  * [42] Jeremy Petch, Shuang Di and Walter Nelson “Opening the Black Box: The Promise and Limitations of Explainable Machine Learning in Cardiology” In _Canadian Journal of Cardiology_ 38, 2022, pp. 204–213
  * [43] James Robins “A New Approach to Causal Inference in Mortality Studies with a Sustained Exposure Period: Application to Control of the Healthy Worker Survivor Effect” In _Mathematical Modelling_ 7, 1986, pp. 1393–1512
  * [44] Erik Sachs, Veronica Jackson and Ulrik Sartipy “Household Disposable Income and Long-Term Survival after Pulmonary Resections for Lung Cancer” In _Thorax_ 75.9, 2020, pp. 764–770
  * [45] Harald Sanftmann “3D Visualization of Multivariate Data”, 2012
  * [46] Aaron L. Sarvet, Kerollos Nashat Wanis, Mats J. Stensrud and Miguel A. Hernán “Graphical Description of Partial Exchangeability” In _Epidemiology_ 31, 2020, pp. 365–368
  * [47] Yu-Ming Shen, Lien D. Le, Rory Wilson and Ulrich Mansmann “Graphical Presentation of Patient-Treatment Interaction Elucidated by Continuous Biomarkers: Current Practice and Scope for Improvement” In _Methods of Information in Medicine_ 56.1, 2017, pp. 13–27
  * [48] Andrew M. Smith, John P. Christodouleas and Wei-Ting Hwang “Understanding the Predictive Value of Continuous Markers for Censored Survival Data Using a Likelihood Ratio Approach” In _BMC Medical Research Methodology_ 19.108, 2019
  * [49] Tyler J. VanderWeele “On Well-Defined Hypothetical Interventions in the Potential Outcomes Framework” In _Epidemiology_ 29.4, 2018, pp. e24–e25
  * [50] Xiaofei Wang, Fangfang Bai, Herbert Pang and Stephen L. George “Bias-Adjusted Kaplan-Meier Survival Curves for Marginal Treatment Effect in Observational Studies” In _Journal of Biopharmaceutical Statistics_ 29.4, 2019, pp. 592–605
  * [51] Meagan M. Wasfy and Aaron L. Baggish “Exercise Dose in Clinical Practice” In _Circulation_ 133.23, 2016, pp. 2297–2313
  * [52] Daniel Westreich and Stephen R. Cole “Positivity in Practice” In _American Journal of Epidemiology_ 171.6, 2010, pp. 674–677
  * [53] Claus O. Wilke “Fundamentals of Data Visualization: A Primer on Making Informative and Compelling Figures” O’Reilly, 2019
  * [54] Hui Yang, Rui Tang, Mike Hale and Jing Huang “A Visualization Method Measuring the Performance of Biomarkers for Guiding Treatment Decisions” In _Pharmaceutical Statistics_ 15, 2016, pp. 152–164
  * [55] Jessica G. Young, Mats J. Stensrud, Eric J. Tchetgen Tchetgen and Miguel A. Hernán “A Causal Framework for Classical Statistical Estimands in Failure Time Settings with Competing Events” In _Statistics in Medicine_ 39.8, 2020, pp. 1199–1236
  * [56] Daniella A. Zipkin, Craig A. Umscheid, Nancy L. Keating, Elizabeth Allen, KoKo Aung, Rebecca Beyth, Scott Kaatz, Devin M. Mann, Jeremy B. Sussman, Deborah Korenstein, Connie Schardt, Avishek Nagi, Richard Sloane and David A. Feldstein “Evidence-Based Risk Communication: A Systematic Review” In _Annals of Internal Medicine_ 161.4, 2014, pp. 270–280

##  Appendix A Alternative Graphics

As discussed in the main text, the three-dimensional survival surface
estimated using g-computation can be visualized in many different ways. The
manuscript focuses on the _survival area plot_ and the _survival contour
plot_. Below we give a short description of the other possibilities and
illustrate them, again using data from the getABI study. Table 2 gives a
concise listing of the most important properties of the different
visualization techniques. Because some of the plots are based on summary
statistics, we first define those quantities.

Using the same notation as in the main text, we can define the _counterfactual
survival time quantile_ as:

$Q_{z}(p)=min\left(t|\hat{S}_{z}(t)\leq p\right),$ (5)

where $p$ is some value between 0 and 1, specifying the quantile of interest.
This quantity can be interpreted as the time until $p\cdot 100$ percent of all
individuals in the target-population would have died if they had received
treatment $Z=z$.

The _counterfactual restricted mean survival time_ on the other hand can be
defined as:

$RMST_{z}(\lambda)=\int_{0}^{\lambda}\hat{S}_{z}(t)dt,$ (6)

where $\lambda$ is some value greater than 0. It can be interpreted as the
average time free from an event until time point $\lambda$ that would be
observed if all individuals in the target population had received treatment
$Z=z$.

Table 2: Properties of different visualization techniques based on estimates of $\hat{S}_{z}(t)$. Plot | $x$-axis | $y$-axis | $z$-axis | Color | Allows | Allows
---|---|---|---|---|---|---
CI | non-linear
| effects
Survival Area | $t$ | $\hat{S}_{z}(t)$ | – | $Z$ | No | No
Survival Contour | $t$ | $Z$ | – | $\hat{S}_{z}(t)$ | No | Yes
Survival Heatmap | $t$ | $Z$ | – | $\hat{S}_{z}(t)$ | No | Yes
3D Survival Surface | $t$ | $\hat{S}_{z}(t)$ | $Z$ | – | No | Yes
Value Specific | $t$ | $\hat{S}_{z}(t)$ | – | Possibly $Z$ | Yes | Yes
Survival Curve
Landmark Survival | $Z$ | $\hat{S}_{z}(t)$ | – | Possibly $t$ | Yes | Yes
Probability
Survival Time Quantile | $Z$ | $\hat{Q}_{z}(p)$ | – |  Possible $p$ | Yes | Yes
Restricted Mean | $Z$ | $\widehat{RMST}_{z}(\lambda)$ | – | Possibly $\lambda$ | Yes | Yes
Survival Time
Figure 4: Survival heatmaps for the causal effect of the ABI on survival,
created using four different Cox proportional hazards regression models and
g-computation. Figure 5: Categorized Survival heatmaps for the causal effect
of the ABI on survival, created using four different Cox proportional hazards
regression models and g-computation. Figure 6: The survival probability at
different points in time $t$ as a function of the ABI, created using four
different Cox proportional hazards regression models and g-computation. Figure
7: The survival probability at $t=6$ with bootstrap confidence intervals as a
function of the ABI, created using four different Cox proportional hazards
regression models and g-computation. The confidence interval was calculated
using the bootstrap percentile method with $1000$ bootstrap samples. Figure 8:
Survival time quantiles as a function of the ABI, created using four different
Cox proportional hazards regression models and g-computation. Figure 9:
Restricted mean survival times as a function of the ABI, created using four
different Cox proportional hazards regression models and g-computation. Figure
10: Survival curves for different values of the ABI, created using four
different Cox proportional hazards regression models and g-computation. Figure
11: Survival 3D surface plot displaying the effect of the ABI on the survival
probability, created using model (2) from table 1 and g-computation. Figure
12: Survival 3D surface plot displaying the effect of the ABI on the survival
probability, created using model (3) from table 1 and g-computation. Figure
13: Survival 3D surface plot displaying the effect of the ABI on the survival
probability, created using model (4) from table 1 and g-computation. Figure
14: Survival 3D surface plot displaying the effect of the ABI on the survival
probability, created using model (5) from table 1 and g-computation.

##  Appendix B Estimating & Depicting Causal Contrasts

Instead of relying on the estimand defined in equation 1 in the main text
directly, it may be useful to depict causal contrasts based on it instead. One
alternative is to use the difference or ratio between the counterfactual
survival probability at $t$ if all individuals had received treatment $Z$, and
the counterfactual survival probability at $t$ if all individuals had received
a specific reference value of $Z$ as treatment, denoted as $\tau$. The
difference is then defined as:

$\Delta_{\tau}(t,z)=S_{\tau}(t)-S_{z}(t),$ (7)

and the ratio is defined as:

$\phi_{\tau}(t,z)=\frac{S_{\tau}(t)}{S_{z}(t)}.$ (8)

In some cases there may be good reasons to choose a specific $\tau$. For
example, if the main goal was to estimate the causal effect of a drug dosage
on the survival probability, using a dosage of $\tau=0$ may be a reasonable
strategy. If that is not the case, one may instead use the survival
probability at $t$ under “treatment-as-usual” conditions. We previously
defined the actually observed survival time as $T$. This probability is then
defined as:

$S(t)=P(T>t).$ (9)

Therefore the causal difference is defined as:

$\Delta_{KM}(t,z)=S(t)-S_{z}(t),$ (10)

and the ratio is defined as:

$\phi_{KM}(t,z)=\frac{S(t)}{S_{z}(t)}.$ (11)

In practice, $S(t)$ can be estimated from the crude sample using standard
techniques, such as the Kaplan-Meier estimator if $T$ is subject to right-
censoring.

Note that these contrast-based target estimands correspond to different causal
questions than the original target estimand defined by equation 1 in the main
text. The causal question associated with the original target estimand can be
translated to be “If every individual in the target population had received
treatment $Z=z$, what proportion of individuals would have survived up to and
including time $t$?”. Table 3 shows the causal question associated with each
of the target estimands.

Table 3: The specific causal questions associated with each of the target estimands. Target Estimand |  Causal Question
---|---
$S_{z}(t)$ |  If every individual in the target population had received treatment $Z=z$, what proportion of individuals would have survived up to and including time $t$?
$\Delta_{\tau}(t,z)$ |  What is the difference between the proportion of individuals that would have survived up to and including time $t$ if every person in the target population had received treatment $Z=\tau$ and the proportion of individuals that would have survived up to and including time $t$ if every person in the target population had received treatment $Z=z$?
$\phi_{\tau}(t,z)$ |  What is the ratio between the proportion of individuals that would have survived up to and including time $t$ if every person in the target population had received treatment $Z=\tau$ and the proportion of individuals that would have survived up to and including time $t$ if every person in the target population had received treatment $Z=z$?
$\Delta_{KM}(t,z)$ |  What is the difference between the proportion of individuals that have survived up to and including time $t$ under the ’treatment-as-usual’ regime and the proportion of individuals that would have survived up to and including time $t$ if every person in the target population had received treatment $Z=z$?
$\phi_{KM}(t,z)$ |  What is the ratio between the proportion of individuals that have survived up to and including time $t$ under the ’treatment-as-usual’ regime and the proportion of individuals that would have survived up to and including time $t$ if every person in the target population had received treatment $Z=z$?

Tables 4–7 displays estimates of these quantities for $Z=0.3$ and $Z=1$ at
$t=4$ and $t=6$, estimated using all of the models described in the main text.
The associated standard errors and 95% confidence intervals were calculated
using the bootstrap percentile method with $1000$ bootstrap samples.

Table 4: Estimates of $\Delta_{KM}$ in the getABI data for some combinations of $t$ and $z$ and all models discussed in the main text. Confidence intervals and standard errors were calculated using $1000$ bootstrap samples. Model | $t$ | $z$ | $\hat{\Delta}_{KM}$ | $SE\left(\hat{\Delta}_{KM}\right)$ |  $CI_{95}$ (Lower) |  $CI_{95}$ (Upper)
---|---|---|---|---|---|---
Linear (unadjusted) | 4 | 0.4 | 0.127 | 0.051 | 0.038 | 0.234
Linear (unadjusted) | 4 | 1.0 | -0.007 | 0.004 | -0.016 | -0.001
Linear (unadjusted) | 6 | 0.4 | 0.210 | 0.076 | 0.069 | 0.355
Linear (unadjusted) | 6 | 1.0 | -0.011 | 0.006 | -0.024 | -0.002
Linear (adjusted) | 4 | 0.4 | 0.088 | 0.041 | 0.017 | 0.173
Linear (adjusted) | 4 | 1.0 | -0.006 | 0.004 | -0.014 | -0.000
Linear (adjusted) | 6 | 0.4 | 0.145 | 0.063 | 0.029 | 0.271
Linear (adjusted) | 6 | 1.0 | -0.009 | 0.005 | -0.021 | -0.001
B-Splines (unadjusted) | 4 | 0.4 | 0.209 | 0.095 | 0.019 | 0.414
B-Splines (unadjusted) | 4 | 1.0 | -0.019 | 0.008 | -0.035 | -0.003
B-Splines (unadjusted) | 6 | 0.4 | 0.329 | 0.127 | 0.036 | 0.562
B-Splines (unadjusted) | 6 | 1.0 | -0.032 | 0.013 | -0.058 | -0.003
B-Splines (adjusted) | 4 | 0.4 | 0.176 | 0.071 | 0.033 | 0.319
B-Splines (adjusted) | 4 | 1.0 | -0.017 | 0.008 | -0.033 | -0.002
B-Splines (adjusted) | 6 | 0.4 | 0.274 | 0.097 | 0.061 | 0.448
B-Splines (adjusted) | 6 | 1.0 | -0.027 | 0.013 | -0.052 | -0.003
Table 5: Estimates of $\phi_{KM}$ in the getABI data for some combinations of $t$ and $z$ and all models discussed in the main text. Confidence intervals and standard errors were calculated using $1000$ bootstrap samples. Model | $t$ | $z$ | $\hat{\phi}_{KM}$ | $SE\left(\hat{\phi}_{KM}\right)$ |  $CI_{95}$ (Lower) |  $CI_{95}$ (Upper)
---|---|---|---|---|---|---
Linear (unadjusted) | 4 | 0.4 | 1.168 | 0.084 | 1.044 | 1.368
Linear (unadjusted) | 4 | 1.0 | 0.992 | 0.004 | 0.982 | 0.999
Linear (unadjusted) | 6 | 0.4 | 1.373 | 0.207 | 1.095 | 1.870
Linear (unadjusted) | 6 | 1.0 | 0.987 | 0.007 | 0.970 | 0.997
Linear (adjusted) | 4 | 0.4 | 1.111 | 0.060 | 1.020 | 1.250
Linear (adjusted) | 4 | 1.0 | 0.993 | 0.004 | 0.984 | 1.000
Linear (adjusted) | 6 | 0.4 | 1.231 | 0.132 | 1.039 | 1.544
Linear (adjusted) | 6 | 1.0 | 0.988 | 0.007 | 0.973 | 0.999
B-Splines (unadjusted) | 4 | 0.4 | 1.312 | 0.215 | 1.022 | 1.895
B-Splines (unadjusted) | 4 | 1.0 | 0.979 | 0.009 | 0.961 | 0.997
B-Splines (unadjusted) | 6 | 0.4 | 1.742 | 0.666 | 1.047 | 3.655
B-Splines (unadjusted) | 6 | 1.0 | 0.961 | 0.016 | 0.930 | 0.996
B-Splines (adjusted) | 4 | 0.4 | 1.250 | 0.133 | 1.039 | 1.564
B-Splines (adjusted) | 4 | 1.0 | 0.981 | 0.009 | 0.964 | 0.997
B-Splines (adjusted) | 6 | 0.4 | 1.547 | 0.329 | 1.084 | 2.377
B-Splines (adjusted) | 6 | 1.0 | 0.966 | 0.015 | 0.937 | 0.996
Table 6: Estimates of $\Delta_{0.6}$ in the getABI data for some combinations of $t$ and $z$ and all models discussed in the main text. Confidence intervals and standard errors were calculated using $1000$ bootstrap samples. Model | $t$ | $z$ | $\hat{\Delta}_{0.6}$ | $SE\left(\hat{\Delta}_{0.6}\right)$ |  $CI_{95}$ (Lower) |  $CI_{95}$ (Upper)
---|---|---|---|---|---|---
Linear (unadjusted) | 4 | 0.4 | 0.055 | 0.025 | 0.014 | 0.109
Linear (unadjusted) | 4 | 1.0 | -0.079 | 0.029 | -0.137 | -0.025
Linear (unadjusted) | 6 | 0.4 | 0.086 | 0.034 | 0.026 | 0.154
Linear (unadjusted) | 6 | 1.0 | -0.134 | 0.047 | -0.223 | -0.046
Linear (adjusted) | 4 | 0.4 | 0.037 | 0.019 | 0.006 | 0.080
Linear (adjusted) | 4 | 1.0 | -0.057 | 0.025 | -0.108 | -0.011
Linear (adjusted) | 6 | 0.4 | 0.058 | 0.027 | 0.011 | 0.115
Linear (adjusted) | 6 | 1.0 | -0.096 | 0.041 | -0.175 | -0.021
B-Splines (unadjusted) | 4 | 0.4 | 0.116 | 0.090 | -0.053 | 0.311
B-Splines (unadjusted) | 4 | 1.0 | -0.112 | 0.037 | -0.190 | -0.041
B-Splines (unadjusted) | 6 | 0.4 | 0.171 | 0.120 | -0.090 | 0.408
B-Splines (unadjusted) | 6 | 1.0 | -0.190 | 0.059 | -0.315 | -0.077
B-Splines (adjusted) | 4 | 0.4 | 0.111 | 0.068 | -0.030 | 0.247
B-Splines (adjusted) | 4 | 1.0 | -0.082 | 0.032 | -0.147 | -0.024
B-Splines (adjusted) | 6 | 0.4 | 0.162 | 0.094 | -0.046 | 0.334
B-Splines (adjusted) | 6 | 1.0 | -0.138 | 0.051 | -0.239 | -0.040
Table 7: Estimates of $\phi_{0.6}$ in the getABI data for some combinations of $t$ and $z$ and all models discussed in the main text. Confidence intervals and standard errors were calculated using $1000$ bootstrap samples. Model | $t$ | $z$ | $\hat{\phi}_{0.6}$ | $SE\left(\hat{\phi}_{0.6}\right)$ |  $CI_{95}$ (Lower) |  $CI_{95}$ (Upper)
---|---|---|---|---|---|---
Linear (unadjusted) | 4 | 0.4 | 1.073 | 0.041 | 1.017 | 1.170
Linear (unadjusted) | 4 | 1.0 | 0.911 | 0.033 | 0.845 | 0.972
Linear (unadjusted) | 6 | 0.4 | 1.153 | 0.091 | 1.036 | 1.369
Linear (unadjusted) | 6 | 1.0 | 0.829 | 0.060 | 0.717 | 0.943
Linear (adjusted) | 4 | 0.4 | 1.047 | 0.028 | 1.007 | 1.116
Linear (adjusted) | 4 | 1.0 | 0.936 | 0.028 | 0.879 | 0.987
Linear (adjusted) | 6 | 0.4 | 1.092 | 0.056 | 1.015 | 1.228
Linear (adjusted) | 6 | 1.0 | 0.877 | 0.052 | 0.777 | 0.974
B-Splines (unadjusted) | 4 | 0.4 | 1.174 | 0.181 | 0.937 | 1.662
B-Splines (unadjusted) | 4 | 1.0 | 0.876 | 0.041 | 0.789 | 0.955
B-Splines (unadjusted) | 6 | 0.4 | 1.384 | 0.493 | 0.881 | 2.836
B-Splines (unadjusted) | 6 | 1.0 | 0.764 | 0.072 | 0.616 | 0.903
B-Splines (adjusted) | 4 | 0.4 | 1.158 | 0.118 | 0.964 | 1.443
B-Splines (adjusted) | 4 | 1.0 | 0.909 | 0.035 | 0.835 | 0.973
B-Splines (adjusted) | 6 | 0.4 | 1.325 | 0.266 | 0.939 | 2.007
B-Splines (adjusted) | 6 | 1.0 | 0.827 | 0.062 | 0.706 | 0.950

It is difficult to use these point estimates without any visualizations,
because they are dependent on both $z$ and $t$. The proposed graphics can be
adapted in a straightforward fashion to display estimates of
$\Delta_{\tau},\Delta_{KM},\phi_{\tau}$ and $\phi_{KM}$ instead of estimates
of $S_{z}(t)$ by simply substituting $\hat{S}_{z}(t)$ with one of those
quantities. Below we give examples for all of these quantities estimated from
the getABI study using some of the previously shown types of graphics.

Figure 15: Continuously and discretely colored survival area plots displaying
the causal effect of the ABI on survival using $\hat{\Delta}_{KM}$, created
using different Cox proportional hazards regression models and g-computation.
A: Using model 1 (unadjusted), continuous color shading; B: Using model 1
(unadjusted), discrete color shading; C: Using model 2 (adjusted), continuous
color shading; D: Using model 2 (adjusted), discrete color shading. Figure 16:
Continuously and discretely colored survival area plots displaying the causal
effect of the ABI on survival using $\hat{\Delta}_{0.6}$, created using
different Cox proportional hazards regression models and g-computation.
A: Using model 1 (unadjusted), continuous color shading; B: Using model 1
(unadjusted), discrete color shading; C: Using model 2 (adjusted), continuous
color shading; D: Using model 2 (adjusted), discrete color shading. Figure 17:
Continuously and discretely colored survival area plots displaying the causal
effect of the ABI on survival using $\hat{\phi}_{KM}$, created using different
Cox proportional hazards regression models and g-computation.
A: Using model 1 (unadjusted), continuous color shading; B: Using model 1
(unadjusted), discrete color shading; C: Using model 2 (adjusted), continuous
color shading; D: Using model 2 (adjusted), discrete color shading. Figure 18:
Continuously and discretely colored survival area plots displaying the causal
effect of the ABI on survival using $\hat{\phi}_{0.6}$, created using
different Cox proportional hazards regression models and g-computation.
A: Using model 1 (unadjusted), continuous color shading; B: Using model 1
(unadjusted), discrete color shading; C: Using model 2 (adjusted), continuous
color shading; D: Using model 2 (adjusted), discrete color shading. Figure 19:
Survival contour plots displaying the causal effect of the ABI on survival
using $\hat{\Delta}_{KM}$, created using four different Cox proportional
hazards regression models and g-computation. Figure 20: Survival contour plots
displaying the causal effect of the ABI on survival using
$\hat{\Delta}_{0.6}$, created using four different Cox proportional hazards
regression models and g-computation. Figure 21: Survival contour plots
displaying the causal effect of the ABI on survival using $\hat{\phi}_{KM}$,
created using four different Cox proportional hazards regression models and
g-computation. Figure 22: Survival contour plots displaying the causal effect
of the ABI on survival using $\hat{\phi}_{0.6}$, created using four different
Cox proportional hazards regression models and g-computation. Figure 23:
$\hat{\Delta}_{KM}$ at different points in time $t$ as a function of the ABI,
created using four different Cox proportional hazards regression models and
g-computation. Figure 24: $\hat{\Delta}_{0.6}$ at different points in time $t$
as a function of the ABI, created using four different Cox proportional
hazards regression models and g-computation. Figure 25: $\hat{\phi}_{KM}$ at
different points in time $t$ as a function of the ABI, created using four
different Cox proportional hazards regression models and g-computation. Figure
26: $\hat{\phi}_{0.6}$ at different points in time $t$ as a function of the
ABI, created using four different Cox proportional hazards regression models
and g-computation.

##  Appendix C Survival Area Plots with Non-Monotonic Effects

As described in the main text, using a standard survival area plot it is
impossible to depict non-monotonic effects of the continuous variable on the
survival probability. Whenever there is a curved relationship between the
continuous variable and the survival probability, some parts of the areas to
be drawn will occupy the same space on the two-dimensional representation,
making the plot unusable. In our illustrative example this is the case in
model (4) and model (5), in which B-splines were used to model the effect of
the ABI.

The only possible way to still use survival area plots in this case is to
divide the plot into multiple facets. Each facet only contains a survival area
plot for different ranges of the covariate space in which its effect on the
survival probability is either non-increasing or non-decreasing. Figure 27
shows the resulting faceted plot when using model (4) and figure 28 shows the
results when using model (5).

Figure 27: Faceted survival area plot displaying the causal effect of the ABI
on survival as estimated using model (4). Figure 28: Faceted survival area
plot displaying the causal effect of the ABI on survival as estimated using
model (5).

The left facet of figure 27 contains the survival area plot for values of the
ABI between $0.35$ and $1.16$, while the right facet contains the survival
area plot for the rest of the considered range of the ABI. This works because
the survival probability strictly increases with higher ABI values in the left
facet and it strictly decreases with higher ABI values in the right facet. If
the relationship between the variable of interest and the survival probability
is more complex, the plot would have to be divided into more facets. Plots
with such relationships become difficult to understand very fast, which is why
we recommend using different plotting strategies, such as survival contour
plots, in those cases instead.

It is important to note that this strategy of dividing the survival area plot
into different segments only works if the shape of the relationship between
the continuous variable and the survival probability is the same at each point
in time. If there are interactions between time and the shape of the
covariate, this strategy is not viable. Additionally, it requires the survival
probability to be strictly non-increasing with increasing time over the whole
range of the continuous variable.

##  Appendix D On the Assumption of Random Right-Censoring

Throughout the article, we explicitly assumed that right-censoring is
independent of any other variable. This assumption allowed us to ignore
censoring in the definition of the target estimand and the description of the
causal identifiability assumptions. In reality, however, multiple variables
may have complex relationships with the occurrence of censoring. In this
section, we give a short description about the implications of dependent
right-censoring on the target estimand, required assumptions and estimation
methodology, based on work by Hernán and Robins (2020).

First, the target estimand needs to be slightly changed in this case. We are
not just interested in the population-level counterfactual survival
probability if the value of $Z$ was set to $z$ for every individual in the
population. Instead we are interested in the population-level counterfactual
survival probability if the value of $Z$ was set to $z$ for every individual
in the population _and_ if they had remained uncensored. This is technically
what we are interested in whether right-censoring is random or not, but if it
is completely random it can be safely ignored. Let $C$ be a dichotomous
indicator of right-censoring with $C=0$ indicating the absence of censoring
and $C=1$ indicating that an observation is censored. Keeping the rest of the
notation of section 2.1, the target estimand can be re-defined as:

$S_{z}(t)=E\left(I\left(T^{Z=z,C=0}>t\right)\right).$ (12)

Essentially, this means that we are no longer interested in the causal effect
of a single treatment on $Z$, but instead that we are interested in the causal
effect of the _joint treatment_ on $Z$ and $C$.

This changes the assumptions required to make $S_{z}(t)$ identifiable. We
still need to make the _no interference_ , the _counterfactual consistency_
assumption, the _conditional exchangeability_ assumption and the _positivity_
assumption, but those are no longer made on just $Z$, but instead have to be
made for the joint treatment on $Z$ and $C$. This makes the estimation of
$S_{z}(t)$ using g-computation more complicated, especially when the
underlying causal structure is highly complex. Westreich et al. (2012)
contains an applied example and more discussion about this topic.

Literature:
Miguel A. Hernán and James M. Robins. _Causal Inference: What If_. CRC Press,
2020. [section 8.4 and 17]
Daniel Westreich, Stephen R. Cole, Jessica G. Young, Frank Palella, Phyllis C.
Tien, Lawrence Kingsley, Stephen J Gange, Miguel A. Hernán. “The parametric
g-formula to estimate the effect of highly active antiretroviral therapy on
incident AIDS or death”. In: _Statistics in Medicine_ 31.18 (2012).

##  Appendix E Additional Results

The following text includes a histogram of the ABI in the used getABI cohort,
as well as the results of graphical tests of the proportional hazards
assumption for all Cox models used and graphics used to identify the correct
functional form of the ABI on survival.

###  E.1 Distribution of the ABI

Figure 29: Histogram showing the distribution of the ABI in the used
population.

###  E.2 Graphical Tests of the Proportional Hazards Assumption

Figure 30: A graphical test for the proportional hazards assumption for each
variable in Cox model (1). The blue lines are non-parametric locally weighted
regressions. Figure 31: A graphical test for the proportional hazards
assumption for the Ankle-Brachial-Index in Cox model (2). The blue line is a
non-parametric locally weighted regression. Figure 32: A graphical test for
the proportional hazards assumption for each variable in Cox model (3). The
blue lines are non-parametric locally weighted regressions. Figure 33: A
graphical test for the proportional hazards assumption for the Ankle-Brachial-
Index in Cox model (4). The blue line is a non-parametric locally weighted
regressions. Figure 34: A graphical test for the proportional hazards
assumption for each variable in Cox model (5). The blue lines are non-
parametric locally weighted regressions.

###  E.3 On the Functional Form of the Relationship Between ABI and Survival

To choose the appropriate functional form of the relationship between the ABI
and the survival time in our applied example, we relied on different graphical
methods. First, we fitted an intercept-only Cox model and plotted the
martingale residuals of this model against the observed values of the ABI
using a LOESS smoother, as shown in figure 35, which clearly indicated a non-
linear relationship. Based on this graphical assessment we tried out different
specifications, using simple polynomial terms (figures 36 and 37), using
b-splines (figures 38–40) and using natural splines (figures 41–44) and
compared the results graphically using survival contour plots. All of those
methods resulted in clear non-linear relationships.

Figure 35: Martingale Residuals of an intercept-only Cox model versus the ABI.
The blue line shown is a non-parametric locally weighted regression plus its
95% confidence band. Figure 36: Survival contour plot created using a Cox
model including the ABI and ABI2. Figure 37: Survival contour plot created
using a Cox model including the ABI, ABI2 and ABI3. Figure 38: Survival
contour plot created using a Cox model including the ABI modelled using
B-splines with three degrees of freedom. Figure 39: Survival contour plot
created using a Cox model including the ABI modelled using B-splines with four
degrees of freedom. Figure 40: Survival contour plot created using a Cox model
including the ABI modelled using B-splines with five degrees of freedom.
Figure 41: Survival contour plot created using a Cox model including the ABI
modelled using natural splines with two degrees of freedom. Figure 42:
Survival contour plot created using a Cox model including the ABI modelled
using natural splines with three degrees of freedom. Figure 43: Survival
contour plot created using a Cox model including the ABI modelled using
natural splines with four degrees of freedom. Figure 44: Survival contour plot
created using a Cox model including the ABI modelled using natural splines
with five degrees of freedom.