# Analysing the radio flux density profile of the M31 galaxy: a possible dark
matter interpretation

Man Ho Chan1 , Chu Fai Yeung1, Lang Cui2 , Chun Sing Leung3
1 Department of Science and Environmental Studies, The Education University of
Hong Kong, Tai Po, Hong Kong, China
2 Xinjiang Astronomical Observatory, Chinese Academy of Sciences, Urumqi,
China
3 Hong Kong Polytechnic University, Hong Kong, China
chanmh<EMAIL_ADDRESS>
(Accepted XXXX, Received XXXX)

###### Abstract

Some recent studies have examined the gamma-ray flux profile of our Galaxy to
determine the signal of dark matter annihilation. However, the results are
controversial and no confirmation is obtained. In this article, we study the
radio flux density profile of the M31 galaxy and show that it could manifest a
possible signal of dark matter annihilation. By comparing the likelihoods
between the archival observed radio flux density profile data and the
predicted radio flux density profile contributed by dark matter and stellar
emission, we can constrain the relevant dark matter parameters. Specifically,
for the thermal annihilation cross section via the $b\bar{b}$ channel, the
best-fit value of dark matter mass is $\sim 30$ GeV, which is consistent with
the results of many recent studies. We expect that this method would become
another useful way to constrain dark matter, which is complementary to the
traditional radio analyses and the other indirect detections.

###### keywords:

(Cosmology:) dark matter; radio continuum: galaxies

††pagerange: Analysing the radio flux density profile of the M31 galaxy: a
possible dark matter interpretation–References††pubyear: XXXX

##  1 Introduction

In the past several decades, observational data of galaxies and galaxy
clusters revealed the existence of dark matter. For example, observations
indicate that most galactic rotation curves are flat in the galactic outskirt
regions, which contradict to the theoretical predictions based on the luminous
mass distributions in galaxies. Also, the dynamical masses of most galaxy
clusters probed from their hot gas distributions are much larger than their
total luminous masses. Therefore, many astrophysicists believe that a large
amount of unknown particles called dark matter exist in the universe. However,
none of the particles in the Standard Model can satisfy the unique properties
of dark matter - almost no interaction with ordinary matter except gravity.

Many theoretical models suggest that dark matter would self-annihilate to give
high-energy gamma rays, electrons, positrons and neutrinos. In particular,
some recent analyses of gamma-ray and radio observations have claimed the
discoveries of potential signals of annihilating dark matter. For example,
using the gamma-ray energy spectrum of our Galaxy or some nearby globular
clusters (e.g. 47 Tuc, Omega Centauri) observed from the Fermi-Large Area
Telescope, some studies show that dark matter with mass $m\approx 30-40$ GeV
annihilating via the $b\bar{b}$ quark channel can best explain the observed
spectral shapes (Daylan et al., 2016; Calore et al., 2015; Abazajian & Keeley,
2016; Cholis, Linden & Hooper, 2019; Brown et al., 2018, 2019). The
constrained annihilation cross sections are surprisingly close to the thermal
annihilation cross section $\langle\sigma v\rangle\approx 2.2\times 10^{-26}$
cm3 s-1 predicted by standard cosmology (Steigman, Dasgupta & Beacom, 2012).

On the other hand, indirect detection of dark matter using radio waveband is
also good for constraining dark matter parameters. For example, some earlier
studies have used the radio data of our Galaxy to constrain dark matter
(Blasi, Olinto & Tyler, 2003; Aloisio, Blasi & Olinto, 2004; Borriello, Cuoco
& Miele, 2009). Later studies have started to focus on the Milky Way satellite
galaxies such as Large Magellanic Cloud (LMC) (Tasitsiomi, Siegal-Gaskins &
Olinto, 2004; Baltz & Wai, 2004; Siffert et al., 2011), nearby galaxies (e.g.
M31 and M33 galaxies) (Borriello et al., 2010; Egorov & Pierpaoli, 2013; Chan,
2018) and galaxy clusters (Colafrancesco, Profumo & Ullio, 2006; Chan et al.,
2020b; Chan & Lee, 2019, 2021) to search for dark matter signals. Some recent
studies show that certain ranges of dark matter mass can give very good fits
to the radio continuum spectra of some galaxy clusters (Chan & Lee, 2019,
2021).

In particular, most of the previous radio studies have used the total
integrated radio emission flux (Blasi, Olinto & Tyler, 2003; Aloisio, Blasi &
Olinto, 2004), radio intensity of certain ideal regions (Borriello et al.,
2010; Siffert et al., 2011) or total integrated radio frequency (energy)
spectrum (multi-frequency approach) (Tasitsiomi, Siegal-Gaskins & Olinto,
2004; Borriello, Cuoco & Miele, 2009; Chan, 2018) of a structure to constrain
dark matter. Nevertheless, using the radial emission profile can also provide
useful constraints for dark matter. For example, some previous studies have
shown that the gamma-ray emissions in our Galactic Center trace the dark
matter density distribution, which provide an indirect potential evidence of
dark matter annihilation (Daylan et al., 2016; Calore et al., 2015). However,
on the contrary, some later studies show that the gamma-ray emissions may
probably trace the stellar density profile rather than the dark matter density
profile (Bartels et al., 2018). Therefore, whether the gamma-ray flux profile
follows the dark matter distribution in our Galactic Center has now become a
controversial issue.

Unfortunately, the resolution of gamma-ray telescopes is not high enough for
us to analyse the central gamma-ray flux profiles of other extragalactic
targets. Therefore, searching dark matter signals by using the gamma-ray flux
profiles of other galaxies is very difficult. Nevertheless, current very large
radio telescopes or radio interferometers can obtain radio data with a very
high resolution and sensitivity. The resolution can be as small as 10” so that
we can observe the radio flux density profiles of galaxies or galaxy clusters.
Therefore, radio flux density profile analysis can be used to serve as an
independent and complementary study to search for dark matter annihilation
signals. In view of this, however, not so much attention has been paid in
considering the radial emission profile (the functional form) in radio bands
to constrain dark matter. In this article, we use the radio flux density
profile of the M31 galaxy as an example to illustrate how this can be done. We
show that some potential dark matter annihilation signals can be manifested in
the radio flux profile of the M31 galaxy. The best-fit value of $m$ for
thermal dark matter generally agrees with the popular range suggested by some
gamma-ray studies.

##  2 Extracting the radio flux profile data

We revisit the archival radio data of the central region of the M31 galaxy
reported in Gießübel & Beck (2014). The data were obtained by the combination
of the Very Large Array (VLA) and the Effelsberg radio telescope data, which
consist of two different frequencies ($\nu=4.85$ GHz and $\nu=8.35$ GHz). The
useful data cover sizes of the central region $r\leq 1.8$ kpc ($\nu=4.85$ GHz)
and $r\leq 1.2$ kpc ($\nu=8.35$ GHz). The resolution for both observing
frequencies is 15”.

We extract the radio flux density per beam size for different positions (RA
and DEC) from the radio map published. We take the coordinates
RA=00h42m44.3503s and DEC=41∘16’09” to be the center of the M31 galaxy (Evans
et al., 2010). As the dark matter distribution and the baryonic mass
distribution in the bulge region are close to spherically symmetric, we take
the azimuthal averaging of the radio flux density (per beam size) in
concentric bins for different angular radii $\theta$ from the center. The
angular distance for each bin is taken as 26” ($1.27\times 10^{-4}$ rad),
which is equivalent to 0.1 kpc at the M31 center (assuming distance to the
galaxy $D=785$ kpc (Egorov & Pierpaoli, 2013)). As there are fluctuations in
the radio flux density for the same angular radius, we take the $1\sigma$
standard deviation as our radio flux density profile uncertainties. Using this
method, we can get the radio flux profiles (radio flux density per beam size
against $\theta$) of the central region of the M31 galaxy for two radio
frequencies (see Fig. 1). The solid angle of the beam size is
$\Delta\Omega=1.66\times 10^{-8}$ sr.

##  3 The dark matter annihilation model

Dark matter annihilation would give a large amount of high-energy electrons
and positrons. These electrons and positrons would emit synchrotron radiation
in radio bands when there is a large magnetic field strength. The power for
synchrotron emission (with energy $E$) takes the following form (Aloisio,
Blasi & Olinto, 2004; Profumo & Ullio, 2010):

$P_{DM}(E,\nu,\vec{r})=\frac{\sqrt{3}e^{3}}{m_{e}c^{2}}B(\vec{r})F(\nu/\nu_{c}),$
(1)

where $\nu_{c}$ is the critical synchrotron frequency, $B$ is the magnetic
field strength and $F$ is the synchrotron kernel function.

Beside synchrotron radiation cooling, the high-energy electrons and positrons
would also cool down via inverse Compton scattering (ICS), Bremsstrahlung
radiation and Coulomb loss (Colafrancesco, Profumo & Ullio, 2006; Egorov &
Pierpaoli, 2013). In particular, synchrotron cooling and the ICS cooling would
dominate at the M31 galactic center. Due to the very high cooling rate, most
of the high-energy electrons and positrons produced from dark matter
annihilation would quickly lose their energy before leaving the central
region. However, the actual diffusion length of the high-energy electrons and
positrons depends on the diffusion models. For the outer part of a galaxy
(e.g. outside the bulge region), the scale of the magnetic field
irregularities is comparable to the gyroradius $r_{g}$ of the electrons and
positrons. As a result, the effect of turbulent diffusion is important so that
the diffusion of the electrons and positrons is quite efficient in the outer
part of a galaxy. Nevertheless, in the central region of a galaxy, the picture
is much more uncertain and the scale of magnetic irregularities would be much
smaller. Therefore, the diffusion coefficient is significantly reduced (Regis
& Ullio, 2008). Like our Galaxy, the diffusion in the deep central region of
the M31 galaxy can be well described by the Bœhm diffusion model, which means
that diffusion effect is not very significant (Regis & Ullio, 2008). The
effective diffusion length of the Bœhm diffusion model can be estimated by the
simple random walk model (Bœhm et al., 2004). The stopping distance of an
electron with energy 1 GeV is $d_{s}\sim\sqrt{r_{g}ct_{c}}\sim 1$ pc, where
$t_{c}\sim 10^{15}$ s is the cooling time. Therefore, almost all of the
electrons and positrons emitted due to dark matter annihilation would be
confined in the central region.

By considering the diffusion-cooling equation, as the diffusion term is not
important, the equilibrium energy spectrum of the high-energy electrons and
positrons become (Borriello, Cuoco & Miele, 2009; Storm et al., 2013; Egorov &
Pierpaoli, 2013):

$\frac{dn_{e}}{dE}=\frac{\langle\sigma
v\rangle\rho_{DM}^{2}}{2m^{2}b(E)}\int_{E}^{m}\frac{dN_{e}}{dE^{\prime}}dE^{\prime},$
(2)

where $\rho_{DM}$ is the density profile of dark matter, $b(E)$ is the cooling
rate and $dN_{e}/dE^{\prime}$ is the injected energy spectrum due to dark
matter annihilation, which can be obtained numerically for different
annihilation channels (Cirelli et al., 2011). The radio emissivity contributed
by dark matter annihilation is mainly determined by the peak frequency so that
the line-of-sight synchrotron radiation flux density (the radio flux density)
within a solid angle $\Delta\Omega$ can be simplified by combining Eq. (1) and
Eq. (2) with the monochromatic approximation (Bertone et al., 2009; Profumo &
Ullio, 2010):

$S_{DM}\approx\frac{1}{4\pi\nu}\left[\frac{9\sqrt{3}\langle\sigma
v\rangle}{2m^{2}(1+C)}E(\nu)Y(\nu,m)\right]\int\rho_{DM}^{2}ds\Delta\Omega,$
(3)

where $s$ is the line-of-sight distance defined as
$r=\sqrt{D^{2}+s^{2}-2Ds\cos\theta}$,
$Y(\nu,m)=\int_{E(\nu)}^{m}(dN_{e}/dE^{\prime})dE^{\prime}$ and
$E(\nu)=14.6(\nu/{\rm GHz})^{1/2}(B/{\rm\mu G})^{-1/2}$ GeV. Here, $C$ is the
correction factor of the ICS cooling contribution. The radiation density due
to synchrotron cooling is $\omega_{B}\approx 9.5$ eV/cm3. For the ICS
radiation density, the value for the M31 galaxy central region has not been
well-determined. Nevertheless, we know that the value for the Milky Way centre
is about $\omega_{i}=8$ eV/cm3 (Egorov & Pierpaoli, 2013). Also, the mass of
the M31 bulge is about 2.3 times of the mass of the Milky Way bulge (Blaña et
al., 2018; Karukes et al., 2020). Therefore, assuming that they have similar
luminosity-to-mass ratio, the ICS radiation density for the M31 galaxy central
region should be $\omega_{i}\approx 19$ eV/cm3. Therefore, we have $C\approx
2$ (i.e. ICS cooling is two times larger than the synchrotron cooling). For
the central magnetic field, previous observations show that the magnetic field
is close to a constant $B=15-19$ $\mu$G for $r=0.2-1$ kpc (Gießübel & Beck,
2014). By taking the mean magnetic field strength $B=17$ $\mu$G, the energies
at peak frequencies are $E(4.85~{}\rm GHz)=7.8$ GeV and $E(8.35~{}\rm
GHz)=10.2$ GeV. Besides, standard cosmology predicts that $\langle\sigma
v\rangle=2.2\times 10^{-26}$ cm3/s for $m>10$ GeV if dark matter is thermally
produced (Steigman, Dasgupta & Beacom, 2012). We will first follow this
standard paradigm to reduce our free parameters in fitting. As a result,
$S_{DM}$ would be a function of $\nu$, $m$ and $\theta$ only.

##  4 Data analysis

As the magnetic field strength is almost constant in the region $r=0.2-1.0$
kpc only (Gießübel & Beck, 2014), we consider the data from angular radius
$\theta=0.0002-0.0012$ rad (i.e. 0.15-0.95 kpc) for our analysis. The dark
matter density is assumed to follow the Navarro-Frenk-White (NFW) density
profile $\rho_{DM}=\rho_{s}r_{s}/r$ (for small $r$) as recent studies show
that it can give a very good fit for the M31 galaxy rotation curve data with
small uncertainties (Sofue, 2015). The density parameters are
$\rho_{s}=(2.23\pm 0.24)\times 10^{-3}M_{\odot}$ pc-3 and $r_{s}=(34.6\pm
2.1)$ kpc (Sofue, 2015). The line-of-sight integral in Eq. (3) is a function
of the angular radius $\theta$. For very small $\theta$, it can be
analytically written as (first-order approximation in $\theta$):

$I(\theta)=\int\rho_{DM}^{2}ds\approx
2\sqrt{\pi}\Gamma(3/2)\left(\frac{\rho_{s}^{2}r_{s}^{2}}{D\theta}\right).$ (4)

We fit the data with our dark matter model for various annihilation channels.
For each frequency $\nu$, the dark matter mass $m$ would be the only free
parameter in the fitting. Generally speaking, the dark matter model can give
good fits for the flux profile data for both frequencies. The goodness of fits
can be determined by the $\chi^{2}$ value, which is defined as

$\chi^{2}=\sum_{i}\left[\frac{S_{m}(\theta_{i})-S_{o}(\theta_{i})}{\sigma_{o}(\theta_{i})}\right]^{2},$
(5)

where $S_{m}(\theta_{i})$, $S_{o}(\theta_{i})$ and $\sigma_{o}(\theta_{i})$
are the data of the modeled flux density, observed flux density and
uncertainties of the observed flux density at different angular radii
$\theta_{i}$. Specifically, we find that the $b\bar{b}$ channel with $m\sim
25$ GeV can give the best fit (i.e. smallest $\chi^{2}$) among all
representative popular annihilation channels ($\chi^{2}=7.1$) (see Table 1 and
Fig. 2). We also separate the fit for each radio frequency. The corresponding
$\chi^{2}$ values are 2.6 and 4.5 for 4.85 GHz and 8.35 GHz respectively.

We also fit the data with the stellar model. We assume that the normal
astrophysical emissions such as pulsar emissions are tracing the stellar
distribution. The stellar distribution in the bulge is commonly modeled by the
Hernquist profile (Sadoun, Mohayaee & Colin, 2014):

$\rho_{b}\propto\frac{1}{r(r+r_{b})^{3}},$ (6)

where $r_{b}=0.61$ kpc for the M31 galaxy. The radio emission due to
conventional astrophysical processes (e.g. pulsar emissions) $S_{*}$ is
proportional to the stellar density profile. Therefore, we define the line-of-
sight integral for the stellar radio emission as
$\tilde{I}(\theta)=\int[r(r+r_{b})^{3}]^{-1}ds$ and we can write
$S_{*}=k\tilde{I}(\theta)$. The proportionality constant $k$ is a free
parameter to fit the observed flux profile for each frequency. Since there are
two different frequencies in the fits, we have totally two free parameters for
fitting. We find that the stellar model can give an overall better fit for the
observed flux profiles ($\chi^{2}=3.2$) compared with the dark matter model
(see Table 1), especially for the profile of $\nu=8.35$ GHz ($\chi^{2}=0.2$).
The radio flux density profiles for $S_{*}$ are shown in Fig. 3.

Although the stellar model can give an overall better fit compared with the
dark matter model, one interesting feature is that the dark matter model gives
a better fit for the 4.85 GHz data while the stellar model gives a better fit
for the 8.35 GHz. In fact, the radio emission due to dark matter annihilation
via the $b\bar{b}$ channel would be more dominant in the low frequency regime
because its energy spectrum has a steeper spectral index. Our results might
suggest that considering both stellar model and dark matter model would be
able to get better fits for both data at 4.85 GHz and 8.35 GHz.

We consider the two-component model $S_{DM}+S_{*}$ to fit the observed flux
profile data. By varying the free parameters $m$ and $k$, we find that a much
better overall fit can be obtained ($\chi^{2}=0.8$) for $m=30$ GeV (see Table
1). The radio flux profiles of the two-component model and the corresponding
components are shown in Fig. 4. Compared with the stellar model, the two-
component model improves the overall fit with a $1.6\sigma$ statistical
preference. Although it does not indicate a very large signal, including the
dark matter annihilation contribution can get an overall better fit. The ratio
of dark matter contribution to stellar contribution is approximately ranging
from 1:1 to 3:1. Therefore, the central radio flux profile data of the M31
galaxy neither simply trace the stellar profile nor the dark matter density
profile, but better trace the stellar plus dark matter density profile.

We have assumed the thermal annihilation cross section $\langle\sigma
v\rangle=2.2\times 10^{-26}$ cm3 s-1 in the above analysis. However, it is
possible that dark matter was created non-thermally and the cross section
would be different from the standard value. Therefore, we release the cross
section as a free parameter so that more general model-independent best-fit
values of cross section as a function of dark matter mass could be obtained.
In Fig. 5, we plot the best-fit average dark matter annihilation cross section
(via $b\bar{b}$ channel) against dark matter mass ($m=20-200$ GeV) based on
the two radio flux density profile datasets (the black line), in which the
best-fit $\chi^{2}$ values are same as that in the thermal dark matter
scenario ($\chi^{2}=0.8$). Such a relation can give us the best-fit dark
matter mass when we know the values of the annihilation cross section based on
different theoretical models or cosmological parameters. Furthermore, if the
annihilation cross section is very large, the radio flux density profile
contributed by dark matter would probably exceed the flux profile data by a
significant amount. Therefore, we can determine the upper limit of the
annihilation cross section as a function of dark matter mass within a
$2\sigma$ margin. In Fig. 5, we also plot the $2\sigma$ upper limit of the
annihilation cross section for $m=20-200$ GeV (the red line).

Figure 1: The radio flux density profile data of the central region of the M31
galaxy. The original radio map is taken from Gießübel & Beck (2014).

Figure 2: The radio flux density profiles of the dark matter annihilation
model (upper: 8.35 GHz; lower: 4.85 GHz). Here $m=25$ GeV for the $b\bar{b}$
channel and $m=185$ GeV for the $e^{+}e^{-}$ channel.

Figure 3: The radio flux density profiles of the stellar model (upper: 8.35
GHz; lower: 4.85 GHz).

Figure 4: The radio flux density profiles of the two-component model (upper:
8.35 GHz; lower: 4.85 GHz). Here we have assumed $m=30$ GeV annihilating via
the $b\bar{b}$ channel.

Figure 5: The best-fit (black line) and the 2$\sigma$ upper limit (red line)
of the annihilation cross section against dark matter mass. Here we have
assumed dark matter annihilating via the $b\bar{b}$ channel.
Table 1: The fitting properties for different models. Model | Channel | $m$ | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}$
---|---|---|---|---|---
|  | (GeV) | (4.85 GHz) | (8.35 GHz) | (both)
$S_{*}$ | - | - | 2.9 | 0.2 | 3.2
$S_{DM}$ | $b\bar{b}$ | 25 | 2.6 | 4.5 | 7.1
| $e^{+}e^{-}$ | 185 | 5.2 | 4.0 | 9.2
$S_{DM}+S_{*}$ | $b\bar{b}$ | 30 | 0.5 | 0.3 | 0.8

##  5 Discussion

In this article, we discuss an analysis of dark matter annihilation by using
the central radio flux density profile of a galaxy. With the observed central
radio flux density data of the M31 galaxy at two frequencies, we show that
adding dark matter contribution to the background radio emission due to
stellar component (the two-component model) can give an overall better fit for
the data. The statistical preference of the two-component model is $1.6\sigma$
compared with the stellar model (the null hypothesis). Therefore, it reveals
some potential signal of dark matter annihilation. However, the data and the
available frequencies are quite limited so that the resulting flux density
profiles have large uncertainties. We expect that a larger set of data with
higher resolution and more observing frequencies can provide a much better
analysis to confirm the annihilation signal.

In our analysis, with the thermal annihilation cross section, the best-fit $m$
is $\sim 30$ GeV for the two-component model (via $b\bar{b}$ channel).
Surprisingly, this value and the annihilation channel is consistent with many
recent studies using gamma-ray (Daylan et al., 2016; Brown et al., 2018, 2019)
and radio data (Chan & Lee, 2021). We have assumed that the annihilation cross
section follows the thermal annihilation cross section predicted by standard
cosmology (Steigman, Dasgupta & Beacom, 2012). Therefore, thermally produced
dark matter with mass $\sim 30$ GeV annihilating via $b\bar{b}$ channel has
become one of the most probable sets of parameters for dark matter. Further
observations and analyses are required to confirm the above claim.
Nevertheless, if dark matter was not created thermally, a large range of best-
fit dark matter mass is still possible to account for the radio flux density
profile of the M31 galaxy.

Previous studies using gamma-ray data of our Galaxy have performed similar
analyses. Some studies claim that the gamma-ray flux profile traces dark
matter distribution while other studies claim that the gamma-ray flux profile
traces stellar distribution. Therefore, it is still a controversial issue.
However, the resolution of current gamma-ray detectors is not high enough to
get the central gamma-ray flux profiles of other nearby galaxies. Fortunately,
the resolution of current radio telescopes is able to fulfill the task. If we
can observe the radio flux density profiles of some nearby galaxies with high
resolution and different frequencies, some better analyses or clearer signals
could be obtained to verify our results. In fact, many previous radio studies
mainly focus on constraining dark matter by using the radio intensity of a
region (e.g. radio intensity of M33 galaxy or LMC, see Borriello et al.
(2010); Siffert et al. (2011)) or radio frequency spectrum (multi-frequency
approach) (e.g. Tasitsiomi, Siegal-Gaskins & Olinto (2004); Chan (2018)) of
galaxies or galaxy clusters. They could identify some ideal regions with the
lowest signal to noise ratio for dark matter detection in the radio wavebands.
Some other studies have obtained the radio sky map of our Galaxy and used the
sky map data to constrain dark matter (Borriello, Cuoco & Miele, 2009).
However, these studies have not explicitly examined the likelihoods between
the observed radial emission profile and the predicted radial emission profile
contributed by dark matter annihilation. In our study, we show that using high
resolution radio density flux profile at the central region of a galaxy is
good for constraining dark matter. Therefore, observing and analysing the
radio flux density profile (i.e. the radial emission profile) would be another
important way to detect the signal of dark matter annihilation and constrain
dark matter properties, which are complementary to the traditional radio
analyses (using total integrated radio flux or frequency spectrum) (Blasi,
Olinto & Tyler, 2003; Aloisio, Blasi & Olinto, 2004; Tasitsiomi, Siegal-
Gaskins & Olinto, 2004; Borriello, Cuoco & Miele, 2009; Chan, 2018), cosmic-
ray (including gamma-ray) analyses (Ackermann et al., 2015; Abdallah et al.,
2016; Ambrosi et al., 2017; Aguilar et al., 2019) and neutrino analyses
(Albert et al., 2020) of dark matter annihilation.

##  6 Acknowledgements

The work described in this paper was supported by a grant from the Research
Grants Council of the Hong Kong Special Administrative Region, China (Project
No. EdUHK 28300518) and the Internal Research Fund from The Education
University of Hong Kong (RG 2/2019-2020R). Lang Cui thanks for the support
from the National Key R&D Program of China (No. 2018YFA0404602), the National
Natural Science Foundation of China (NSFC grant No. 61931002 & U1731103) and
the Youth Innovation Promotion Association of the CAS (No. 2017084).

##  7 Data availability statement

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abazajian & Keeley (2016) Abazajian K. N. & Keeley R. E., 2016, Phys. Rev. D 93, 083514.
  * Abdallah et al. (2016) Abdallah H. et al. [H.E.S.S. Collaboration], 2016, Phys. Rev. Lett. 117, 111301.
  * Ackermann et al. (2015) Ackermann M. et al., 2015, Phys. Rev. Lett. 115, 231301.
  * Aguilar et al. (2019) Aguilar M. et al., 2019, Phys. Rev. Lett. 122, 041102.
  * Albert et al. (2020) Albert A. et al., 2020, arXiv:2003.06614.
  * Aloisio, Blasi & Olinto (2004) Aloisio R., Blasi P. & Olinto A. V., 2004, J. Cosmol. Astropart. Phys. 05, 007.
  * Ambrosi et al. (2017) Ambrosi G. et al., 2017, Nature 552, 63.
  * Baltz & Wai (2004) Baltz E. A. & Wai L., 2004, Phys. Rev. D 70, 023512.
  * Bartels et al. (2018) Bartels R., Storm E., Weniger C. & Calore F., 2018, Nature Astron. 2, 819.
  * Bertone et al. (2009) Bertone G., Cirelli M., Strumia A. & Taoso M., 2009, J. Cosmol. Astropart. Phys. 03, 009.
  * Blaña et al. (2018) Blaña Diaz M. et al., 2018, Mon. Not. R. Astron. Soc. 481, 3210.
  * Blasi, Olinto & Tyler (2003) Blasi P., Olinto A. V. & Tyler C., 2003, Astropart. Phys. 18, 649.
  * Bœhm et al. (2004) Bœhm C., Hooper D., Silk J., Casse M. & Paul J., 2004, Phys. Rev. Lett. 92, 101301.
  * Borriello, Cuoco & Miele (2009) Borriello E., Cuoco A. & Miele G., 2009, Phys. Rev. D 79, 023518.
  * Borriello et al. (2010) Borriello E., Longo G., Miele G., Siffert B. B., Tabatabaei F. S. & Beck R., 2010, Astrophys. J. 709, L32.
  * Brown et al. (2018) Brown A. M., Lacroix T., Lloyd S., Bœhm C. & Chadwick P., 2018, Phys. Rev. D 98, 041301.
  * Brown et al. (2019) Brown A. M., Massey R., Lacroix T., Strigari L. E., Fattahi A. & Bœhm C., arXiv:1907.08564.
  * Calore et al. (2015) Calore F., Cholis I., McCabe C. & Weniger C., 2015, Phys. Rev. D 91, 063003.
  * Chan (2018) Chan M. H., 2018, Mon. Not. R. Astron. Soc. 474, 2576.
  * Chan & Lee (2019) Chan M. H. & Lee C. M., 2019, Phys. Dark Univ. 26, 100355.
  * Chan & Lee (2021) Chan M. H. & Lee C. M., 2021, Mon. Not. R. Astron. Soc. 500, 5583 (arXiv:1912.03640).
  * Chan et al. (2020b) Chan M. H., Lee C. M., Ng C.-Y. & Leung C. S., 2020b, Astrophys. J., 900, 126 (arXiv:2007.06547).
  * Cholis, Linden & Hooper (2019) Cholis I., Linden T. & D. Hooper, 2019, Phys. Rev. D 99, 103026.
  * Cirelli et al. (2011) Cirelli M. et al., 2011, J. Cosmol. Astropart. Phys. 03, 051.
  * Colafrancesco, Profumo & Ullio (2006) Colafrancesco S., Profumo S. & Ullio P., 2006, Astron. Astrophys. 455, 21.
  * Colafrancesco, Marchegiani & Beck (2015) Colafrancesco S., Marchegiani P. & Beck G., 2015, J. Cosmol. Astrophys. Phys. 02, 032.
  * Daylan et al. (2016) Daylan T., Finkbeiner D. P., Hooper D., Linden T., Portillo S. K. N., Rodd N. L. & Slatyer T. R., 2016, Phys. Dark Univ. 12, 1.
  * Egorov & Pierpaoli (2013) Egorov A. E. & Pierpaolo E., 2013, Phys. Rev. D 88, 023504.
  * Evans et al. (2010) Evans I. N. et al., 2010, Astrophys. J. Supp. 189, 37.
  * Gießübel & Beck (2014) Gießübel R. & Beck R., 2014, Astron. Astrophys. 571, A61.
  * Karukes et al. (2020) Karukes E. V., Benito M., Iocco F., Trotta R. & Geringer-Sameth A., 2020, J. Cosmol. Astropart. Phys. 05, 033.
  * Profumo & Ullio (2010) Profumo S. & Ullio P., Particle Dark Matter: Observations, Models and Searches, ed. G. Bertone, Cambridge: Cambridge University Press, chapter 27 (2010).
  * Regis & Ullio (2008) Regis M. & Ullio P., 2008, Phys. Rev. D 78, 043505.
  * Sadoun, Mohayaee & Colin (2014) Sadoun R., Mohayaee R. & Colin J., 2014, Mon. Not. R. Astro. Soc. 442, 160.
  * Siffert et al. (2011) Siffert B. B., Limone A., Borriello E., Longo G. & Miele G., 2011, Mon. Not. R. Astron. Soc. 410, 2463.
  * Sofue (2015) Sofue Y., 2015, Publ. Astron. Soc. Jpn. 67, 759.
  * Steigman, Dasgupta & Beacom (2012) Steigman G., Dasgupta B. & Beacom J. F., 2012, Phys. Rev. D 86, 023506.
  * Storm et al. (2013) Storm E., Jeltema T. E., Profumo S. & Rudnick L., 2013, Astrophys. J. 768, 106.
  * Tasitsiomi, Siegal-Gaskins & Olinto (2004) Tasitsiomi A., Siegal-Gaskins J. M. & Olinto A. V., 2004, Astropart. Phys. 21, 637.


secttoc2

# The nph2ph-transform: applications to the statistical analysis of completed
clinical trials

Sean M. Devlin1 and John O’Quigley2

(1Memorial Sloan Kettering Cancer Center, New York, USA
2Department of Statistical Science, University College London, U.K. )

###### Abstract

We present several illustrations from completed clinical trials on a
statistical approach that allows us to gain useful insights regarding the time
dependency of treatment effects. Our approach leans on a simple proposition:
all non-proportional hazards (NPH) models are equivalent to a proportional
hazards model. The nph2ph transform brings an NPH model into a PH form. We
often find very simple approximations for this transform, enabling us to
analyze complex NPH observations as though they had arisen under proportional
hazards. Many techniques become available to us, and we use these to
understand treatment effects better.

Brownian motion, clinical trials, model fit, proportional hazards, time-
dependent effects

##  1 Introduction

There are two main goals to a clinical research study that involves survival
as a main endpoint. The first goal, an important one although relatively easy
to express, is to decide if some experimental treatment results in prolonging
survival. The log-rank test plays a key role in such studies and typically
forms the basis for sample size calculations given a Type 1 error and hoped-
for power under an alternative hypothesis that the experimental treatment does
have an impact. To carry out such calculations, some simplifying assumptions
need to be made, the most common being that if we can reject the null
hypothesis, the alternative will be well approximated by a proportional
hazards model. The great majority of studies make such an assumption and will
sometimes include in the protocol some formal way of testing the validity of
the assumption.

The second main goal will be addressed once the trial and the main results are
completed. The subsequent analysis is no longer restricted by the design
structure, the hypotheses of interest, and the desire to maintain careful
control on statistical error. We are at liberty to carry out a broader
investigation of the treatment effect. While conclusions from such secondary
investigations will lack the impact concerning the main objectives and
endpoints, they can still be of considerable value in model building. The
purpose of model building may be to gain greater insight and understanding of
the myriad of mechanisms that can affect survival or to obtain some pointers
into ways in which we might refine and improve our initial design. This can
help plan future related studies.

Our focus here is on this second goal. There is a rich literature on
prognostic studies and model building (Harrell et al., 1996, Gerds et al.,
2008, Steyerberg et al., 2013) and this can be an important component of the
retrospective analysis of a clinical trial. We do not consider this here since
our focus is somewhat narrower. We consider the study itself, most often
involving two comparative groups, and investigate the impact of the
proportional hazards assumption and to what extent predictive power can be
increased when survival is more accurately modeled by a non-proportional
hazards assumption.

Figure 1: Kaplan-Meier survival curves from a randomized trial in melanoma
(Long et al., 2017). The figure on the top shows the estimated conditional
survival functions obtained under the best-fitting proportional hazards model.
The bottom figure shows the estimated conditional survival functions under the
best-fitting piecewise proportional hazards model after implementing a nph2ph
transform.

A recent study (Figure 1) published in the New England Journal of Medicine
showed the clear advantage of dabrafenib combined with trametinib over placebo
in melanoma patients with a BRAF mutation. The outcome was relapse-free
survival, and, in this regard, the superiority of the treatment was clear. An
assumption of PH looks to be a very reasonable one. However, the PH model is
not a good fit, as shown later in Section 7.1. The simplest possible
approximation, as discussed in Section 4, to the unknown nph2ph transform - a
piecewise PH model with a single changepoint - dramatically improves fit. This
is immediately visible in Figure 1 where not only does the transformed model
provide a much more accurate summary of the Kaplan-Meier curves, but we also
note that there is minimal room for further improvement.

The strength of the treatment effect, measured under an incorrect PH
assumption, will underestimate the true strength of the effect when treatment
effects are time-dependent (Chauvel and O’Quigley, 2017). The bias can be
large, and we will see many situations in which the actual treatment effect
can be double, or even triple, what we had thought it to be due to a poor
model specification. This is a very important finding in its own right and
enables us to provide a more relevant summary of the trial’s outcome. The fact
that a trial’s design may have led to the underestimation of treatment effect,
along with an indication of how much, can be a valuable consideration in
deciding how best to move forward. We mention in passing that the power of
many statistical tests, the log-rank test, in particular, depends on the PH
assumption. Under a more accurate model specification, we may have achieved
greater power.

The nph2ph transform is very much a data analytic tool. By this, we mean that
we put less emphasis on theoretical population models and the resulting
distribution theory and more on the data. Our starting point is the observed
data rather than some postulated model. We introduce an index of fit, which
measures the degree to which observed data can be taken to have a PH form. Our
interest is not in testing the fit - although this could be done - but in
deciding how close to or far from a PH structure the data are. And, in using
the nph2ph transformation, we can make NPH observations appear more like PH
observations.

##  2 Basic methodology

###  2.1 Proportional and non-proportional hazards models

The random variables of interest are the failure time, $T_{i}$, the censoring
time, $C_{i},$ and the treatment indicator covariate, $Z_{i}(\cdot)$,
$i=1,...,n.$ Generally, the treatment indicator, $Z_{i}(t),$ will remain
constant in time. Allowing for time dependency creates no additional
methodological difficulties; in certain situations, for example, cross-over
trials, this may be used to advantage. We use $\mathscr{I}(A)$ as an indicator
function so that, $\mathscr{I}(A)=1$ when $A$ is true and is zero otherwise.
Let $F(t)=P(T<t)$. For each subject $i$ we observe $X_{i}=\min(T_{i},C_{i})$,
and $\delta_{i}=\mathscr{I}(T_{i}\leq C_{i})$. The “at risk” indicator
$Y_{i}(t)$ is defined as, $Y_{i}(t)=\mathscr{I}(X_{i}\geq t).$ The counting
process $N_{i}(t)$ is defined as, $N_{i}(t)=\mathbf{1}\\{T_{i}\leq t,T_{i}\leq
C_{i}\\}.$ $\bar{N}(t)=\sum_{1}^{n}N_{i}(t)$ and we transform the failure
times to the interval (0,1) via the homomorphism $\phi(t),$ described below.
Closely related to $Z$ is the observation at $t$, ${\cal
Z}_{t}=\sum_{i=1}^{n}Z_{i}\left(X_{i}\right)\mathscr{I}\\{\displaystyle\phi_{n}(X_{i})=t,\delta_{i}=1\\},\quad
t\in[0,1],$ in words, an outcome that assumes the value of the covariate
(group assignment) associated with the failure, when there is one, at time
$t,$; otherwise, it assumes the value zero. O’Quigley (2003) and Chauvel and
O’Quigley (2014) chose to suppress the random nature of $T,$ treating it as an
indexing variable $(1/k_{n},2/k_{n},...)$ on the interval (0,1) so that we can
now view ${\cal Z}_{t}$ as a stochastic process, i.e., a collection of random
variables indexed by $t.$ To complete the characterization of the process, we
associate with each $t,$ the distribution $\pi_{i}(\beta(t),t).$ While any
well-defined distribution could be used, we focus on the specific definition

$\displaystyle\pi_{i}(\beta(t),t)=\frac{Y_{i}(t)\exp\\{\beta(t)Z_{i}(t)\\}}{\sum_{j=1}^{n}Y_{j}(t)\exp\\{\beta(t)Z_{j}(t)\\}},\text{
}i=1,\ldots,n,$ (2.1)

which arises naturally in the setting of proportional and non-proportional
hazards models.

We group the models under as general a heading as possible. The most general
model is then the non proportional hazards model written as

$\lambda(t|Z(t))=\lambda_{0}(t)\exp\\{\beta(t)Z(t)\\},$ (2.2)

where $\lambda(t|\cdot)$ is the conditional intensity function,
$\lambda_{0}(t)$ the baseline hazard and $\beta(t)$ the time-varying
regression effect. The above model becomes a proportional hazards model under
the restriction that $\beta(t)=\beta$, a constant i.e.
$\lambda(t|Z(t))=\lambda_{0}(t)\exp\\{\beta Z(t)\\}.$ In the proportional
hazards model, a strictly monotonically increasing transformation in time does
not affect inference of the regression coefficient $\beta$. Indeed, the time
of death and censored values are not used for inference; only their order of
occurrence counts. The transformation of Chauvel and O’Quigley (2014) will
also conserve the order of observed times while simultaneously leading to
great tractability. By applying the inverse transformation, we can go back to
the initial time scale. More precisely, the transformation of Chauvel and
O’Quigley (2014) can be seen as an isomorphism $\phi_{n}$ on the vector
$X_{1},\dots,X_{n},$ taking us from the observed time outcomes to the interval
(0,1). The observed failures map to the values $(1/k_{n},2/k_{n},...,)$ where
$k_{n}$ equates to the observed number of failures in the simplest case. More
generally, we will ignore failures that contain no information relating to
treatment effect, i.e., the set of the larger failures taking place when the
risk set only contains individuals from the same treatment group.

###  2.2 Treatment effect process

A key ingredient in building the nph2ph transform is the treatment effect
process, which is no more than a special case of the regression effect process
(Chauvel and O’Quigley, 2014, 2017, Flandre and O’Quigley, 2019, O’Quigley,
2021) applied to the situation of two treatment groups. Expectations and
variances with respect to the family of probabilities
$\\{\pi_{i}(\beta(t),t),i=1,\dots,n\\}$ can now be defined, and these are all
that are needed to define the treatment effect process.

###### Definition 1

The moments of $Z$ with respect to the family of probabilities
$\\{\pi_{i}(\beta(t),t),i=1,\dots,n\\}$ and in which $t\in[0,1],$ are given
by,

${\cal E}_{\beta(t)}(Z^{m}\mid
t)=\sum_{i=1}^{n}Z^{m}_{i}\left(\phi_{n}^{-1}(t)\right)\pi_{i}(\beta(t),t),$
(2.3)

providing us with the mean ${\cal E}_{\beta(t)}(Z\mid t)$ and variance, ${\cal
V}_{\beta(t)}(Z\mid t)={\cal E}_{\beta(t)}(Z^{2}\mid t)-{\cal
E}^{2}_{\beta(t)}(Z\mid t).$

We use this mean and variance to construct the following standardized process
that we will refer to as the treatment effect process.

###### Definition 2

Let $j=0,1,\dots,k_{n}$ and $t_{j}=j/k_{n}.$ Given $\beta(t_{j}),$ the
treatment effect process at time $t_{j}$ is written:

$U^{*}_{n}\left(\beta(t_{j}),t_{j}\right)=\dfrac{1}{\sqrt{k_{n}}}\sum_{i=1}^{j}{\cal
V}_{\beta(t_{i})}(Z\mid t_{i})^{-1/2}\left\\{\mathscr{Z}(t_{i})-{\cal
E}_{\beta(t_{i})}(Z\mid t_{i})\right\\}.$ (2.4)

We extend this to the whole $[0,1]$ interval by linear interpolation of the
$k_{n}+1$ random variables $\left\\{U_{n}^{*}(\beta(t_{j}),t_{j}),\
j=0,1,\dots,k_{n}\right\\},$ one of which is degenerate, taking on the value
zero with probability one when $t_{0}=0.$ We then have: the standardized score
process, $U_{n}^{*}(\beta(t),t),\ t\in[0,1]$, evaluated at $\beta(t)$ where
for $j=0,\dots,k_{n}$ and $t\in[t_{j},t_{j+1}]$, is defined on (0,1) by:

$\displaystyle
U^{*}_{n}\left(\beta(t),t\right)=U^{*}_{n}\left(\beta(t_{j}),t_{j}\right)+\left(tk_{n}-j\right)\left\\{U^{*}_{n}\left(\beta(t_{j+1}),t_{j+1}\right)-U^{*}_{n}\left(\beta(t_{j}),t_{j}\right)\right\\}.$
(2.5)

Before establishing some needed theoretical results, we can gain some
intuition from Figure 2 that describes a few situations. The first is the
absence of effect, corresponding to Brownian motion. The other situations
describe proportional and non-proportional hazards scenarios. The following
theorem makes these observations more precise.

###### Theorem 1

Let $\beta_{0}$ be a constant function on the class of possible regression
functions $\mathbb{B}$. Under the non-proportional hazards model with
parameter $\beta(t)$, consider

$A_{n}(t)=\dfrac{1}{k_{n}}{\displaystyle\int_{0}^{t}}\left\\{\beta(s)-\beta_{0}\right\\}\dfrac{\mathscr{V}_{\tilde{\beta}(s)}(Z\mid
s)}{\mathscr{V}_{\beta_{0}}(Z\mid s)^{1/2}}d\bar{N}^{*}(s),\,\quad 0\leq t\leq
1,$ (2.6)

where $\tilde{\beta}$ is in the ball with center $\beta$ and radius
$\sup_{t\in[0,1]}\left|\beta(t)-\beta_{0}\right|$ and $N^{*}(s)$ is the
counting process on the transformed time scale. Then, there exist two
constants $C_{1}(\beta,\beta_{0})$ and $C_{2}\in\mathbb{R}^{*+}$ such that
$C_{1}(\beta,\beta)=1$ and

$U^{*}_{n}(\beta_{0},\cdot)-\sqrt{k_{n}}A_{n}\quad\overset{\cal
L}{\underset{n\rightarrow\infty}{\longrightarrow}}\quad
C_{1}(\beta,\beta_{0})\mathscr{W}(t),$ (2.7)

where $\mathscr{W}(t)$ is a standard Brownian motion. Furthermore,

$\sup_{0\leq t\leq
1}\left|A_{n}(t)-C_{2}\int_{0}^{t}\left\\{\beta(s)-\beta_{0}\right\\}ds\right|\quad\overset{P}{\underset{n\rightarrow\infty}{\longrightarrow}}\quad
0.$ (2.8)

###### Corollary 1

Under the non-proportional hazards model with parameter $\beta(t)$,

$U^{*}_{n}\left(\beta(\cdot),\cdot\right)\quad\overset{\cal
L}{\underset{n\rightarrow+\infty}{\longrightarrow}}\quad\mathscr{W}(t).$

A proof of the theorem is given in Chauvel and O’Quigley (2014). The theorem
implies that the regression effect process will trace out a path that closely
parallels that of $\int_{0}^{t}\beta(u)du.$ Simulations back this up, and it
is straightforward to obtain curves such as those illustrated in Figure 2. On
the original time scale, however, things will not look so clear cut so that,
even in the case of proportional hazards, the process can appear to be non-
linear. That said, on the finite set of time points, the transformation
$\phi(t)$ is an isomorphism so when done with any calculation, we can always
go back and forth between the original and the transformed time scales by
employing $\phi(t)$ and $\phi^{-1}(\cdot).$

Figure 2: The treatment effect process $U^{*}_{n}$ under different scenarios.
The top left corresponds to no treatment effect, and the top right corresponds
to a constant proportional hazards treatment effect. The bottom left and right
correspond to a delayed treatment effect that occurs around the 40th
percentile on the transformed time scale and an early treatment effect up
until the 50th percentile.

###  2.3 The nph2ph-transform

In practical situations, as opposed to theoretical study or simulations, the
nph2ph transform will not be known. However, the treatment effect process,
which is observed, provides direct information on the form of the
transformation. As a consequence of some key results concerning the treatment
effect process, we can construct a simple approximation to the unknown
transform and assess to what degree this empirical transform has worked. We
can be assisted in this by both visual and numerical quantification and clear
guidelines on when to decide that our empirically derived transform is, in
some sense, good enough. The following simple lemma is at the core of our
development and states that every non-proportional hazard model,
$\lambda(t|Z(t))=\lambda_{0}(t)\exp\\{\beta(t)Z(t)\\}$, is equivalent to a
proportional hazards model.

###### Lemma 1

For given $\beta(t)$ and covariate $Z(t)$ there exists a constant $\beta_{0}$
and time dependent covariate $Z^{*}(t)$ so that
$\lambda(t|Z(t))=\lambda_{0}(t)\exp\\{\beta(t)Z(t)\\}=\lambda_{0}(t)\exp\\{\beta_{0}Z^{*}(t)\\}\,.$

There is no loss of generality if we take $\beta_{0}=1.$ The result is
immediate upon introducing $\beta_{0}\neq 0$ and defining a time-dependent
covariate $Z^{*}(t)\equiv Z(t)\beta(t)\beta_{0}^{-1}.$ The important thing to
note is that we have the same $\lambda_{0}(t)$ on either side of the equation
and that, whatever the value of $\lambda(t|Z(t))$, for all values of $t$,
these values are exactly reproduced by either expression, i.e., we have
equivalence. Simple though the lemma is, it has strong and important
ramifications, and we continue to develop them here.

##  3 Making improvements to the model

The most common path to model improvement is by including regressors or
predictor variables that improve the precision of our knowledge of the outcome
variable. That path is well understood and not our focus here. Here, we are
only concerned with ways to improve the model once we have already settled on
the set of regressors or predictor variables. Since generalizations are quite
straightforward, we only deal with the simplest case: a single binary
regressor or treatment variable coded (0,1). Model improvements are then
concentrated on the time effects. These improvements come under two closely
related headings: goodness of fit and predictive strength.

###  3.1 Improving model goodness-of-fit

There is a large body of literature focused on the problem of goodness-of-fit
for proportional hazards models beginning with the test of Cox (1972). Many of
these might be potential candidates for our purpose, although we will only
focus on the approach of O’Quigley (2003). That approach is based on the
regression effect process, the treatment effect process in the case of a
single binary variable, and so is fully in tune with the methods we use in
this work. The regression effect itself, characterized by $\beta(t),$ need not
be constant, which is also something we will need, which would rule out many
other methods. Finally, our purpose is usually not to carry out any formal
test with accurate control of error rates but to gain a visual impression of
how well our model can be taken to accurately summarize time dependency.

The visual summary we are looking for is based on the following. Under the
model, the process $U^{*}_{n}\left(\beta(t),t\right)$ can be treated as though
it were a realization of a Brownian motion. Slutsky’s theorem allows us to
conclude that the same will hold for $U^{*}_{n}\left(\hat{\beta}(t),t\right)$
as long as $\hat{\beta}(t)$ provides a consistent estimate of $\beta(t).$
Realistically, we will not be able to achieve that for general $\beta(t)$ but,
for special cases, for instance, $\beta(t)$ constant or piecewise constant
over a small division of the interval $(0,1),$ this can be manageable. From
this we are able to derive the bridged process
$U^{0}_{n}\left(\beta(t),t\right)$;

###### Definition 3

The bridged process is defined by the transformation

$U^{0}_{n}(\beta(t),t)=U^{*}_{n}(\beta(t),t)-t\times U^{*}_{n}(\beta(t),1)$

###### Lemma 2

The process $U^{0}_{n}(\beta(t),t)$ converges in distribution to the Brownian
bridge, in particular $E\,U^{0}_{n}\left(\beta(t),t\right)=0$ and $\mbox{\rm
Cov}\,\\{U^{0}_{n}(\beta(s),s)$, $U^{0}_{n}(\beta(t),t)\\}$ $=s(1-t).$

Once again, these results continue to hold for $U^{0}_{n}(\hat{\beta}(t),t)$
whenever $\hat{\beta}(t)$ is consistent for $\beta(t).$ The Brownian bridge is
also called tied-down Brownian motion for the obvious reason that at $t=0$ and
$t=1$ the process takes the value 0. So when our model is “correct”, i.e., it
can be viewed as an accurate representation of how the observations are
generated, then this bridged process will start and end at 0 and will show no
obvious drift, appearing to walk randomly around the origin axis. How far it
may plausibly wander can be quantified by using:

$\displaystyle\mbox{\rm Pr}\,\left\\{\sup_{u}|U^{0}_{n}(\hat{\beta}(t),u)|\geq
a\right\\}\approx 2\exp(-2a^{2}),$ (3.9)

which follows as a large sample result from Kolmogorov’s theorem;

$\mbox{\rm Pr}\,\left\\{\sup_{u}|U^{0}_{n}(\hat{\beta}(t),u)|\leq
a\right\\}\to 1-2\sum_{k=1}^{\infty}(-1)^{k+1}\exp(-2k^{2}a^{2})\,,\ \ a\geq
0\,.$

This is an alternating sign series, and therefore, if we stop the series at
$k=2$ then the error is bounded by $2\exp(-8a^{2})$ which for most values of
$a$ that we will be interested in will be small enough to ignore.

Given the value of the process is 0 at $t=0$ and $t=1$, the variance of
$U^{0}_{n}(\hat{\beta}(t),t)$ clearly depends on $t$. As an alternative, a
standardized distribution, defined as
$U^{0}_{n}(\hat{\beta}(t),t)/\sqrt{t(1-t)}$, could be utilized (O’Quigley,
2003). This standardized distribution is not well defined at $t=0$ and $t=1$,
so the interval (0,1) needs to be reduced to $(\epsilon_{1},\epsilon_{2})$,
where $\epsilon_{1}>0$ and $\epsilon_{2}<1$. With these components, we can use
the following approximation by Miller and Siegmund (1982),

$\displaystyle\mbox{\rm
Pr}\,\left\\{\sup_{u}\left|\frac{U^{0}_{n}(\hat{\beta}(t),u)}{\sqrt{u(1-u)}}\right|\leq
a\right\\}\approx\frac{4\phi(a)}{a}+\phi(a)\left(a-\frac{1}{a}\right)\log\left\\{\frac{\epsilon_{2}(1-\epsilon_{1})}{\epsilon_{1}(1-\epsilon_{2})}\right\\},$
(3.10)

where $a\geq 0$ and $\phi(x)$ is the standard normal distribution. With this
result, we can use confidence intervals for the bridge process that go to 0
towards $t=0$ and $t=1$. When the bridged process strays into zones with small
probability under the large sample Brownian bridge assumption, this indicates
that the current model is not adequate. This will also be observable from the
treatment effect process itself and will simultaneously suggest a way for
improvement. This can be seen in the examples below.

###  3.2 Improving predictive strength of the model

The predictive ability of a model is typically quantified by some $R^{2}$
measure of explained variation. The reason for this follows as a consequence
of Chebyshev’s inequality showing that the smaller the variance, the more
clustered, on average, the observations will be around the predicted
conditional mean. The conditional mean is all the more accurately located as
the goodness-of-fit improves. And so the predictive strength of the model and
the goodness-of-fit, while quantifying different things, are closely related.
We explore this further below, in particular, using an important result of
Chauvel and O’Quigley (2017) showing that the best fitting model will be the
closest to the model that generates the observations and that will maximize
the population equivalent of our $R^{2}$ measure.

We will measure predictive strength in two ways. We use $R^{2}$ that indicates
the reduction in explained variance due to the treatment indicator variable
and, as a result of the Chebyshev inequality, the degree by which our
predictability has improved. Secondly, we use $\kappa,$ the concordance
coefficient, that estimates how high the probability is that a random subject
from the treatment group will live longer than one from the reference group.

Note that, for any random variables $Z$ and $T$ having second moments, we have
that

$\text{Var}(Z)=E(\text{Var}(Z|T=t))+\text{Var}(E(Z|T=t)).$

From this the definition of explained variation, $\Omega^{2},$ is immediate
and we write,
$\Omega^{2}(\beta(t))=1-\text{Var}(E(Z|T=t))/\text{Var}(Z)=E\text{Var}((Z|T=t))/\text{Var}(Z).$
In this expression we can replace the expectation operator $E$ by ${\cal E}$
and the variance operator Var by ${\cal V}$ and, since these depend upon
$\pi_{i}(\beta(t),t)$ which in turn depends upon $\beta(t)$ we make this
clear, by writing $\Omega^{2}(\beta(t)),$ emphasizing the dependence on the
model as characterized by $\beta(t).$ We also have:

###### Lemma 3

A consistent estimator of $\Omega^{2}(\beta(t))$ is $R^{2}(\hat{\beta}(t))$ in
which $\hat{\beta}(t)$ is any consistent estimator of $\beta(t),$ in
particular $\hat{\beta},$ the usual partial likelihood estimator under
proportional hazards.

These results were shown in O’Quigley and Flandre (1994), where several
properties were also outlined. These include large sample independence of
$R^{2}$ to an independent censoring mechanism, invariance to linear
transforms, and monotonically increasing transforms on time, i.e., rank
invariance. For our purposes, we need more, and the following results support
our development. Suppose that $\mathscr{B}$ defines a class of possible
regression functions $\beta(t).$ This class can be as large as we wish,
restricted in any way we wish, i.e., monotonic increasing or decreasing, and
not necessarily continuous. It does, though, have to be finite. Then:

###### Theorem 2

Under the non–proportional hazards model (2.2) of parameter $\beta(t)$, we
have

$|R^{2}(\beta(t))-R^{2}(\hat{\beta}(t))|\overset{a.s.}{\underset{n\rightarrow+\infty}{\longrightarrow}}0,$

and with probability one, $\underset{{b(t)\in\mathscr{B}}}{\arg\max}\
\lim_{n\rightarrow+\infty}R^{2}(b(t))=\beta(t).$

The proof is given in Chauvel and O’Quigley (2017). If $\beta(t)$ is the true
regression coefficient, the maximum of $R^{2}$ is then well approximated by
$R^{2}(\hat{\beta}(t))$ for a large enough sample size. Using the theorem
results, the treatment effect process can indicate the form of the regression
effect in relation to time. Tools such as smoothing or the projection on a
basis of functions or kernel estimation allow for other approaches (Cai and
Sun, 2003, Hastie and Tibshirani, 1990, Scheike and Martinussen, 2004). We
have not compared the different tools. As an illustration, Figure 2 shows a
constant effect until time $\tau$ followed by a null effect. This is easily
detectable, especially with moderate and larger sample sizes.

We can assume that $\beta(t)$ can be written $\beta_{0}\times b(t)$ where
$b(t)$ is a member of $\mathscr{B}.$ Given $b(t)$ we replace $\beta(t)$ by
$\hat{\beta}(t)=\hat{\beta}_{0}\times b(t)$ and $\hat{\beta}_{0}$ is obtained
in any of the usual ways, eg via maximization of the partial likelihood (Cox,
1972). The function $b(t)$ can be determined graphically from the treatment
effect process, noting that this function is only specified up to an arbitrary
constant. In a single changepoint model, all that is needed to obtain $b(t)$
is the ratio of the two slopes. More formally, let
$\mathscr{B}=\\{b_{1}(t),\dots,b_{m}(t)\\}$ be a set of $m$ functions from
$[0,1]$ to $\mathbb{R}$. Note that while $m$ is not unbounded, it can be as
large as we wish. The selected regression function $b^{*}(t)$ is such that

$b^{*}(t)=\arg\max_{b(t)\in\mathscr{B}}R^{2}\left(b(t)\right).$

The following theorem gives an equivalence between this maximization problem
when $n\rightarrow\infty$ and a problem of minimization of $L^{2}$ norm.

###### Theorem 3

Under the non–proportional hazards model with regression parameter $\beta(t)$
not necessarily in $\cal B$, asymptotically,
$\beta^{*}(t)=\underset{\alpha(t)\in\mathscr{B}}{\arg\max}\
\underset{n\rightarrow\infty}{\lim}R^{2}\left(\alpha(t)\right)$ is the
solution of

$\beta^{*}(t)=\arg\min_{\alpha(t)\in\mathscr{B}}\left\|\beta(t)-\alpha(t)\right\|_{2}=\arg\min_{\alpha(t)\in\mathscr{B}}\left(\int_{0}^{1}\left(\beta(t)-\alpha(t)\right)^{2}dt\right)^{1/2}.$

Other approaches may be possible but we have not explored them in view of a
lack of results regarding their properties. Several estimators of the
explained variation have been proposed in the literature (Choodari-Oskooei et
al., 2012). Here we have avoided presenting an exhaustive review of these
estimators. Instead, we have focused on the $R^{2}$ coefficient introduced by
O’Quigley and Flandre (1994) since it is built with the same ingredients as
the treatment effect process. Our proposal in this paper leans on several
well-established results.

##  4 Working approximations for the nph2ph-transform

Our goal is not to find the best possible estimate of the nph2ph transform, a
goal that would be very challenging both conceptually and operationally.
Instead, we aim to find a simple approximation to the nph2ph transform that
appears to do well enough. Our definition of “well enough” equates to a
transform that results in a bridged process that does not stray outside the
boundaries established according to Equation (3.10), i.e., in Figure 4 (D) for
an upcoming example, one that lies fully inside the unshaded zone. According
to Chauvel and O’Quigley (2017), once such a situation prevails,, there is
very little room for increasing the model’s predictability unless, of course,
we add in new or further covariates. Briefly, if some simple approximation to
nph2ph brings a bridged process, some of which has strayed into the shaded
areas, into one that is fully contained within the unshaded area, then we can
claim an improvement of the proportional hazards approximation and anticipate
a resulting improvement in predictive power. This has been the case with the
examples shown here, and with many others, we have looked at.

While the treatment effect process conveys the needed information on the shape
of the nph2ph transform, it also contains an amount of irregular noise that we
may wish to eliminate. Drawing a well-fitting curve over the process would
provide a working solution but would not be readily reproducible. Most likely,
there are many possible solutions provided by any established curve fitting
technique. In the next two sections, we describe two of these, both of which
we exploit in our examples. The first solution provides a smooth, continuously
differentiable curve. The second is a rough and ready, albeit powerful,
approximation based on a change-point model.

###  4.1 Legendre polynomials with bounded second derivatives

Recall that we aim to find something that is “good enough” for the data at
hand, meaning that we want some working summary of what has been observed.
Orthogonal polynomials can be a useful way to approximate functions. Their
accuracy will depend on the number of included terms, and, in our situation,
they have some appeal in that such functions are both continuous and
differentiable. They are no more than working approximations, and although
many orthogonal bases arise as solutions to differential equations of
mathematical physics, any attempt to make physical connections might not seem
realistic. That said, orthogonal polynomials will allow us to build our
approximations incrementally so that the higher the order, the better the job
they do. There is, though, a downside to that: the problem of overfitting,
whereby a great deal of flexibility can result in random events and background
noise having too strong an influence on the fit. We can readily control this
by not allowing the fitted curve to turn too quickly, i.e., by putting bounds
on the second derivative. One candidate class of orthogonal polynomials that
comes to mind is the class, $P_{0}(t),$ $P_{1}(t),...,$ known as the Legendre
polynomials. We have that $P_{0}(t)=1$ and $P_{1}(t)=t$ so that these provide
the basis for the absence of regression effect and a proportional hazards
effect. Higher orders obtained from

$(n+1)\times P_{n+1}(t)=(2n+1)\times tP_{n}(t)-n\times P_{n-1}(t)\,,$

bring steadily improving fits. In our limited experience, it will generally
suffice to employ no more than lower-order Legendre polynomials. One helpful
fact is that the derivatives, $\Delta(P_{n})$, of Legendre polynomials can be
expressed in terms of polynomials of lower order. This enables us to obtain a
visual impression of the form of $\beta(t).$ In Figure 4, we see a Legendre
approximation and a piecewise PH approximation, described just below. Both
offer major improvements over the simple PH model.

###  4.2 Piecewise proportional hazards models

A piecewise PH model with a single changepoint would be the first step in
approximating the nph2ph transform. Indeed, in the examples we show and many
others we have looked at, this very simple first step away from PH can be very
effective. Given real observations, there are many ways in which such a model
could be fit. Among these would be the changepoint survival models described
by O’Quigley and Pessione (1991), Contal and O’Quigley (1999), and O’Quigley
and Natarajan (2004). We could also focus our attention on the treatment
effect process itself, which, under a simple assumption of a two-stage PH
model, would provide two regression slopes to be estimated. The technique of
Chow (1960) is a natural contender for this estimation although it would not
guarantee continuity of the slopes at the changepoint and, given that our goal
is to model Brownian motion with drift, this would seem like a weakness. Our
preferred approach is to use Euler’s theorem for the drift parameter of
Brownian motion with linear drift. This provides us with two regression slopes
that are continuous everywhere although, of course, differentiable nowhere in
a large sample sense, as indeed is Brownian motion itself. We will appeal to
known theoretical results of Brownian motion and the Brownian bridge.

Let us take the changepoint ${\tau}$ to be the point at which the slope
changes and our estimate of this is based on optimizing fit so that;
$\tau=\arg\max_{i=1}^{k_{n}}L(s)$ where;

$L(s)=\sum_{i}\left\\{\log\pi_{i}(\beta_{1}(X_{i}),X_{i})I(X_{i}<s)+\log\pi_{i}(\beta_{2}(X_{i}),X_{i})I(X_{i}\geq
s)\right\\}.$

Now that we have our estimate for $\tau,$ we can write down the piecewise
linear function $U^{P}_{n}(\beta(t),t)$ that best fits the observations. Note
that by best fit, we mean the function consequent upon maximizing the
likelihood $L(s)$ across all possible changepoints. Given any changepoint, we
only need to estimate the drift coefficient of Brownian motion for before and
after $\tau.$ We do this by appealing to Euler’s theorem whereby a simple
linear connection of the first point on the process to the last (on either
side of $\tau$) would correspond to the maximum likelihood estimate of the
drift parameter for Brownian motion with linear drift. We then have;

$U^{P}_{n}(\beta(t),t)=\frac{t\times
U^{*}_{n}(\beta(t),\tau)\mathscr{I}(t<\tau)}{\tau}+\frac{[(t-\tau)\times
U^{*}_{n}(\beta(t),1)+U^{*}_{n}(\beta(t),\tau)(1-t)]\mathscr{I}(t\geq\tau)}{1-\tau}$

We can see that this function is piecewise linear and, for our purposes, we
need take note of two slopes, the one before and the one after the changepoint
$\tau.$ This provides us with the form for $\beta(t)$ which for computational
purposes we will re-express in terms of $Z(t).$ Consider then,

###### Lemma 4

The piecewise proportional hazards model with time-dependent coefficient
$\beta(t)$ and changepoint $\tau$ on the transformed scale is equivalent to a
proportional hazards model with time-dependent treatment indicator, $Z^{P}(t)$
in which;

$\displaystyle Z^{P}(t)=\mathscr{I}(t\leq\tau)\times
Z(t)+\frac{\tau\mathscr{I}(t>\tau)}{1-\tau}\left\\{\frac{U^{*}_{n}(\beta(t),1)-U^{*}_{n}(\beta(t),\tau)}{U^{*}_{n}(\beta(t),\tau)}\right\\}\times
Z(t)$ (4.11)

###### Proposition 1

Suppose the observations are generated by the NPH model with coefficient
function $\beta(t).$ The set ${\cal B}_{1}$ is the collection of all
proportional hazards models and the set ${\cal B}_{2},$ the collection of all
piecewise PH models with a changepoint at $\tau.$ Then
$ER^{2}(\beta_{2}(t))\geq ER^{2}(\beta_{1})$ where:

$\beta_{j}(t)=\arg\min_{\alpha(t)\in\mathscr{B}_{j}}\left(\int_{0}^{1}\left(\beta(t)-\alpha(t)\right)^{2}dt\right)^{1/2}\,,\
\ \ j=1,2.$

The proposition is simply a corollary to Theorem 3 and is immediate. It tells
us that our predictive strength, as measured by $R^{2},$ will not decrease in
expectation when we replace the simple proportional hazards model by the
single changepoint piecewise proportional hazards model. Again, this assumes
that we do not change the covariates in the model and how they are identified.
We may, of course, purely for computation, re-express a time-dependent effect
by a time-dependent covariate as in Equation 4.11. No optimality properties
are claimed for the transform, and its success is to be judged by its impact
on fit and the improvement in predictive strength.

Equation 4.11 looks a little involved at first glance but is no more than a
time-dependent rescaling of the covariate based on the ratios of the two
slopes of the process. If these slopes are the same, then the covariate is
unchanged. Note also that any PH model remains a PH model under a linear
transformation of the covariate. The product of the regression coefficient,
$\beta$, and the covariate, i.e., the log relative risk, does not change, and
so, in Equation 4.11, the covariate value before $\tau$ can be taken to be the
original coding, $Z(t).$ We can use this observation to find the piecewise
approximation to the nph2ph transform.

##  5 Coefficient of concordance

For our purposes, we make use of $R^{2}$ in the light of its many desirable
properties, in particular, the properties described in the above subsections,
allowing us to build better-fitting models from poorer-fitting ones. Many
other coefficients, similar in spirit to $R^{2}$ have been suggested in the
literature, see for example Harrell et al. (1996), Gönen and Heller (2005),
Uno et al. (2011), Choodari-Oskooei et al. (2012), and these might also serve
the purpose of our choice here. One logical choice since it directly
translates predictability is the coefficient of concordance, $\kappa(A,B).$
Recall that, for $R^{2}$ and two treatment groups, $A$ and $B$ with survival
times $T_{A}$ and $T_{B},$ $R^{2}$ will tend to the value zero when the group
rankings are random, unrelated to the group status, and will tend to the value
one when the sets of rankings become almost perfectly separated. Moving away
from independence, i.e., $R^{2}=0,$ we know that, in expectation, $R^{2}$
increases as the degree of concordance increases (Xu and O’Quigley, 1999) and
we measure concordance via;

$\kappa(A,B)=1-\kappa(B,A)=\mathrm{Pr}\,(T_{A}>T_{B})$

###### Lemma 5

Suppose that $T_{A}$ relates to $T_{B}$ via a proportional hazards model with
parameter $\beta,$ then:

$\kappa(A,B)=\frac{1}{1+\exp(-\beta)}$

When $T_{A}$ and $T_{B}$ are related via the simple exponential model, we can
obtain the result by the convolution of $T_{A}-T_{B}$ and simple integration.
The result then immediately extends to more general proportional hazards
models. A useful measure of predictive strength is the relative risk
$\psi(A,B)$ where $\psi(A,B)=\kappa(A,B)/(1-\kappa(A,B)).$ When
$\kappa(A,B)=0.5,$ there is no regression effect and $\psi(A,B)=1.$ Suppose
that $\psi(A,B)=1.5.$ We could then say that under treatment, $A,$ a patient
has a 50% greater odds than one under treatment $B$ of being the longer
survivor of the two. The most complex model we consider here is a model with a
single changepoint, $\tau.$

###### Lemma 6

For a piecewise PH model and parameters $\beta_{1}(T<\tau)$ and
$\beta_{2}(T>\tau)$ then;

$\displaystyle\kappa(A,B)=\frac{\\{1-\exp(-\tau)\\}\\{1-\exp(-\tau\beta_{1})\\}}{1+\exp(-\beta_{1})}+\frac{\exp\\{-\tau(1+\beta_{2})\\}}{1+\exp(-\beta_{2})}+\exp(-\tau)\\{1-\exp(-\tau\beta_{1})\\}$

Proof: A simple partition of the outcome space leads to;

$\displaystyle\kappa(A,B)=$ $\displaystyle=$
$\displaystyle\Pr(T_{A}>T_{B}|T_{A}<\tau,T_{B}<\tau)\times\Pr(T_{A}<\tau,T_{B}<\tau)+$
$\displaystyle\Pr(T_{A}>T_{B}|T_{A}>\tau,T_{B}<\tau)\times\Pr(T_{A}>\tau,T_{B}<\tau)+$
$\displaystyle\Pr(T_{A}>T_{B}|T_{A}<\tau,T_{B}>\tau)\times\Pr(T_{A}<\tau,T_{B}>\tau)+$
$\displaystyle\Pr(T_{A}>T_{B}|T_{A}>\tau,T_{B}>\tau)\times\Pr(T_{A}>\tau,T_{B}>\tau).$

Using Lemma 5, standard results of the exponential distribution in conjunction
with the rank invariance of proportional hazards models, the independence of
survival times, the simplification that
$\Pr(T_{A}>T_{B}|T_{A}>\tau,T_{B}<\tau)=1$ and
$\Pr(T_{A}>T_{B}|T_{A}<\tau,T_{B}>\tau)=0$, and gathering together all of the
terms, we have the result. While we can see that the details can quickly
become fastidious, it is clear that the extension to several change-points
raises no conceptual difficulties.

##  6 Using the approximated nph2ph transform to estimate the conditional
survival functions

It is common to estimate the survival functions via non-parametric methods
such as Kaplan-Meier or Nelson-Aalen since parametric models, although
relatively easy to fit, will rarely fit well to real data. A proportional
hazards model in the case of say two groups allows us to obtain estimates that
are similar in spirit but that share information. The groups are not
independent. The estimated regression coefficient of the model will fully
capture the dependency structure of the two groups. In the case of
proportional hazards the two survival functions will look like power
transforms of one another. This relationship breaks down for non-proportional
hazards situations, and the PH fitted models may do a poor job in
approximating the survival curves. This is well illustrated by Figure 1 of the
Introduction. The nph2ph transform is a very effective tool that facilitates
the estimation of the survival function in NPH situations.

The properties of the Nelson-Aalen estimates are asymptotically equivalent to
the Kaplan-Meier estimates. In practical situations, the two estimates tend to
all but coincide, and, in our view, the Nelson-Aalen estimate is a better
choice for the development since the relationship to PH and NPH models is more
transparent. To see this, for two groups, $G=0$ and $G=1,$ the survival
functions given by $S_{G}(t)=\exp\\{-\Lambda(t:G)\\}\,,G=0,1,$ in which,

$\displaystyle\Lambda(t:G)=\,\sum_{i=1}^{n}\frac{\mathscr{I}({X_{i}\leq
t})\delta_{i}\exp\\{\beta(X_{i})G\\}}{\sum_{j=1}^{n}Y_{j}(X_{i})\exp\\{\beta(X_{i})Z_{j}(X_{i})\\}},\text{
}$ (6.12)

When we have a single group, or when $\beta(t)=0,$, then this expression
coincides with the usual Nelson-Aalen estimate, which, in turn, all but
coincides with the Kaplan-Meier estimate. In practice we replace $\beta(t)$ by
$\hat{\beta}(t).$ Not precisely knowing $\beta(t)$ implies that we do not know
precisely nph2ph. However, the most elementary approximation to nph2ph can
work very well. We saw this in the example of the Introduction. We have seen
it in other examples, but the lack of space prevents us from showing it here.
For illustration, we consider one further example in which it is immediately
clear that PH provides a poor model, while the most elementary of nph2ph
approximations, a single change point model, indicates a significantly
improved fit. The improvement, as expressed via the conditional survival
function in Equation (6.12), is so strong that any further improvements will
necessarily be very small (see Figure 3). This example is further discussed in
the following section.

Figure 3: Kaplan-Meier curves from a survival study in colorectal cancer
(André et al., 2020). The top figure shows the estimated conditional survival
functions obtained under the best-fitting proportional hazards model. The
bottom figure shows the estimated conditional survival functions under the
best-fitting proportional hazards model after implementing an nph2ph
transform.

##  7 Recent examples from the literature

The nph2ph transform is demonstrated using multiple published phase III trials
in oncology. The trials were selected to illustrate different types of
treatment effects observed across various clinical settings. As the individual
patient-level data were not directly available for each trial, the published
survival curves were digitized to recreate individual-level data for this
illustration (Guyot et al., 2012, Satagopan et al., 2017).

In each randomized trial, a proportional hazards regression model estimated
the hazard ratio associated with each treatment effect. The goodness-of-fit
was assessed for each model by examining the bridged treatment effect process
outlined in Section 3.1. Confidence bounds were added to each bridge plot
using Equation (3.10) to aid model assessment. Two sets of bounds were added
to each figure (see, for example, plot (B) in Figure 4). If the bridged
process were to stay within the center white area, it would indicate the
current model provides a reasonable fit. The border of the light gray area
represents a 90% confidence bound for the bridge, and the dark gray area
represents a 99.9% confidence bound; if the bridged process were to enter
these shaded areas, it would represent different degrees of model inadequacy.

A Shiny app of the nph2ph transform is available to confirm these results or
run on additional data sets (https://nph2ph.shinyapps.io/nph2ph/).

###  7.1 Dabrafenib and Trametinib Compared to Placebo in Patients with
Melanoma Cancer

Long et al. (2017) conducted a phase III randomized trial to compare
combination therapy with dabrafenib and trametinib to a placebo in melanoma
patients with a BRAF mutation. The trial enrolled 870 patients who were
randomized 1:1 to the two treatment arms. The primary endpoint of the study
was relapse-free survival. At the time of trial reporting, a total of 414
patients died or relapsed.

The estimated Kaplan-Meier relapse-free survival curves are shown in Figure 4
(A) for the digitized data set. Using proportional hazards regression, the
estimated hazard ratio was 0.47 (95% CI: 0.38-0.57). The $R^{2}$ was modest at
0.133.

Figure 4:  (A) Estimated relapse-free survival for patients treated with
dabrafenib and trametinib compared to placebo in Long et al. (2017). (B) The
bridged treatment effect process indicating a major model deviation. (C) The
treatment effect process with the identified changepoint $\tau=0.46$ and the
estimated slopes before and after $\tau$. (D) The bridged process after the
implementation of the changepoint. (E) and (F) are the estimated Legendre
polynomials for the treatment effect process and the subsequent bridged
process.

Investigation of the bridged regression process in Figure 4 (B) indicates a
major deviation from proportional hazards. The maximum of the bridged process
was 4.45, which, using the bounds outlined in Equation (3.10), indicates a
very poor model fit.

Using the treatment effect process provided in Figure 4 (C), a changepoint was
identified at transformed time $t=0.46$, which corresponds to approximately
8.5 months on the original time scale. The ratio of the two estimated slopes
shown in the figure before and after the changepoint is -0.024. This results
in an estimated effect of
$\hat{\beta}(t)=-2.08\times(I\\{t<0.46\\}-0.024\times I\\{t\geq 0.46\\})$.
Interpreting these results, there is a large estimated effect on the rate of
relapse or death when treated with dabrafenib and trametinib compared to
placebo in the first eight months following randomization with an estimated
hazard ratio of 0.12. However, after eight months, the estimated effect
largely disappears with a hazard ratio of 1.02.

The $R^{2}$ for this new model increased considerably to 0.306. The bridged
process for this new model is shown in Figure 4 (D). The bridged process
remained within the denoted bounds; therefore, there was little indication of
model inadequacy after implementation of the single changepoint. With this new
model well approximated by proportional hazards, the concordance coefficient
was estimated to be 0.668.

Despite having no indication of model misspecification after implementing a
single changepoint, a two-changepoint model could be explored, given the
severity of proportional hazards deviation in the original data. This richer
model resulted in a modest improvement to fit, with an improvement in $R^{2}$
from 0.306 to 0.327. Collectively, this indicates that a single changepoint
model on its own was enough to overcome the large model misspecification that
was observed when the original data were fit using proportional hazards
regression.

An alternative approach to approximating the nph2ph transform could be based
on the Legendre polynomials as outlined in Section 4.1. One such polynomial is
provided in Figure 4 (E). The bridged process in the next panel indicates a
reasonable model fit. The resulting $R^{2}$ for this polynomial is 0.317, a
modest improvement over the single changepoint model. The estimated effect
under this approximation is
$\hat{\beta}(t)=1.50\Delta(P_{1}(t))-3.06\Delta(P_{2}(t))+2.62\Delta(P_{3}(t))-0.81\Delta(P_{4}(t))$.

###  7.2 Pembrolizumab Compared to Chemotherapy in Patients with Colorectal
Cancer with High Microsatellite Instability

André et al. (2020) conducted a phase III randomized study of pembrolizumab
compared to chemotherapy among advanced or metastatistc colorectal cancer
patients who have tumors with high microsatellite instability or that are
mismatch repair deficient. A total of 153 patients were randomized to receive
pembrolizumab and 154 were randomized to receive chemotherapy. The primary
endpoint was progression-free survival; the study had a median follow-up of
32.4 months.

Figure 5(A) shows the Kaplan-Meier survival curves for the two treatment
groups. From these curves, there appears to be a deviation from the
proportional hazards, a deviation we can also observe from Figure 5(B) and
(C). The maximum of the bridge was 3.51. In the digitized data set, the hazard
ratio was estimated to be 0.59 (95% CI: 0.44-0.79) and the $R^{2}$ was 0.053.

André et al. (2020) recognized the non-proportionality of the data and
supplemented the proportional hazard regression results with an analysis of
the restricted mean survival time. We include this example to illustrate how
the nph2ph transform can be used as a straightforward approach to interpret
the treatment effect in this setting.

A changepoint of 0.55 was identified using the treatment effect process in
Figure 5(C). On the original time scale, the changepoint corresponds to 4.2
months. Using the ratio of -4.40 for the two estimated slopes, the estimated
treatment effect after implementing the changepoint was
$\hat{\beta}(t)=0.41\times(I\\{t<0.55\\}-4.40\times I\\{t\geq 0.55\\})$, which
shows the benefit of pembrolizumab after the changepoint. The $R^{2}$ of the
new model increased considerably to 0.265.

The improvement in fit is also illustrated by examining best-fitting
proportional hazards models before and after nph2ph implementation, as shown
in Figure 3. A clear improvement is observed.

Similar to the previous section, Legendre polynomials were explored as an
alternative nph2ph transform. This is shown in Figure 5(E), and the bridged
process is shown in Figure 5(F). There appears to be no evidence of model
inadequacy. The estimated effect under this approximation is
$\hat{\beta}(t)=-0.98\Delta(P_{1}(t))+1.00\Delta(P_{2}(t))-1.32\Delta(P_{3}(t))+0.38\Delta(P_{4}(t))$.
The $R^{2}$ for this approach was 0.274, a value close to 0.265 based on the
changepoint model.

Figure 5:  (A) Estimated progression-free survival for colorectal cancer
patients treated with pembrolizumab compared to chemotherapy in André et al.
(2020). (B) A major proportional hazards model deviation was identified from
the bridged process. (C) The treatment effect process with the identified
changepoint $\tau=0.55$ and the estimated slopes before and after $\tau$. (D)
The bridged process after the implementation of the changepoint. (E) and (F)
are the estimated Legendre polynomials for the treatment effect process and
the subsequent bridged process.

###  7.3 Cetuximab Compared to Best Supportive Care Patients with Colorectal
Cancer

A phase III randomized study compared cetuximab to best supportive care for
patients with advanced colorectal cancer expressing EGFR (Jonker et al.,
2007). A total of 287 patients received cetuximab, while 285 received best
supportive care. The investigators found that cetuximab significantly improved
both overall survival and progression-free survival.

The progression-free survival curves for the two treatment arms are shown in
Figure 6 (A) for the digitized data set. Visually, from the Kaplan-Meier
curves, there is a clear delay in the treatment effect occurring around two
months following randomization. The estimated hazard ratio for the treatment
effect was 0.69 (95% CI: 0.58-0.82) using proportional hazards regression. The
$R^{2}$ was 0.03.

Figure 6:  (A) Estimated progression-free survival for patients treated with
cetuximab compared to best supportive care in Jonker et al. (2007). (B) The
bridged treatment effect process indicates a deviation from proportion hazards
with a maximum value of 2.63. (C) The treatment effect process with the
identified changepoint $\tau=0.57$. (D) The bridged process after the
implementation of the changepoint was within the denoted boundaries. (E) and
(F) are the estimated Legendre polynomials for the treatment effect process
and the subsequent bridged process.

The bridged process, shown in Figure 6 (B) indicates deviation from
proportional hazards with a maximum of 2.63. Using the regression effect
process in (C), the nph2ph transform identified 0.57 as the changepoint on the
transformed time scale, corresponding to 1.92 months on the original time
scale. The ratio of the two estimated slopes before and after the changepoint
is 5.24. This results in an estimated treatment effect on the transformed time
scale of $\hat{\beta}(t)=-0.089\times(I\\{t<0.57\\}+5.24\times I\\{t\geq
0.57\\})$. The hazard ratios of 0.91 and 0.63 before and after the changepoint
at 1.9 months from the nph2ph transform align with what is visually observed
in the Kaplan-Meier curves.

The $R^{2}$ for the new model is 0.051. The bridged process in Figure 6 (D)
indicates no further model deviation after implementing the single
changepoint. The concordance is estimated to be 0.558 for the new model.

The alternative polynomial approach to nph2ph and the corresponding bridged
process are shown in Figure 6 (E)-(F). The estimated effect for this
approximation is
$\hat{\beta}(t)=-3.54\Delta(P_{1}(t))+3.17\Delta(P_{2}(t))-2.01\Delta(P_{3}(t))+0.58\Delta(P_{4}(t))(t)$.
The $R^{2}$ for this estimate is 0.077. Interestingly, increasing the order of
the Legendre polynomial to $P_{8}(t)$ only results in an $R^{2}$ improvement
to 0.096 despite aligning more closely with the treatment effect process. This
supports the idea that once the bridged process is within the bounds in Figure
6 (F), there will be only a marginal $R^{2}$ gain in further developing a
richer model.

###  7.4 Ipilimumab Compared to Placebo in Patients with Melanoma

In the final example, Eggermont and colleagues conducted a phase III
randomized trial of ipilimumab compared to placebo in patients with stage III
melanoma (Eggermont et al., 2016). A total of 951 were randomized 1:1 to the
two arms. Over a five-year median follow-up of 5.3 years, the authors found
that ipilimumab significantly improved overall and recurrence-free survival.

The overall survival curves are shown in Figure 7 (A) for the digitized data
set. The estimated hazard ratio for the treatment effect was 0.71 (95% CI:
0.58 - 0.88). The bridged process in (B) and the regression effect process in
(C) collectively indicate no notable deviation from proportional hazards.
Therefore, an nph2ph transform is not required for this example.

Figure 7:  (A) Estimated overall survival for patients treated with ipilimumab
compared to placebo in Eggermont et al. (2016). The bridge (B) and the
treatment effect process (C) indicate no model inadequacy.

##  8 Discussion

Our approach in this work is very much a data analytic one. While our focus
and development are built around the nph2ph transform, we limit our efforts to
estimating this transform, which would, in practice, mean finding a good
estimate of a time-transformed regression function, $\beta(t).$ A lot of work
on finding good estimates of $\beta(t)$ has already been done (Murphy and Sen,
1991, Murphy, 1993) and, rather than build on that, we set a more modest goal.
This goal is not to identify as best we can the transform itself but to find a
ready working approximation to it that we consider to be good enough in
practice. For this to work, we need to define what we mean by “good enough”
and we have taken it to be one that results in a good fit for the model.
Although there are many formal tests of fit available and new tests could be
built out of the bridged transform of the treatment effect process, we believe
that our goal is reasonably well accomplished when this bridged process does
not stray outside the white zone of the figures. Our conclusion on fit would
be based on the “eyeball factor”, does it look good, i.e., does it behave as
we would wish so that the bridged process does not find itself in any shaded
area.

The examples and many others that we have studied tend to suggest that for
situations in which the regression function $\beta(t)$ is either constant or
slowly changing and not too volatile, a simple changepoint model will do
enough to produce an improvement in fit that is good enough. This is backed up
by the theoretical findings in O’Quigley (2003), Chauvel and O’Quigley (2014),
and Chauvel and O’Quigley (2017) indicating that once we have achieved a good
fit, any further improvements in fit will not correspond to noticeable changes
in predictability. Let’s not forget that all of this is conditional upon the
regressors in the model. Moving away from a single binary covariate to more
covariates may, of course, greatly impact predictability. That is another
issue and not one we have looked into here.

In a practical setting, the figures shown in the examples can add much to the
discussion. Poor fit and the potential for significant improvements in
predictability are immediately apparent. This is unlike many current
discussions around non-proportionality that take the observed survival
estimates, either the Kaplan-Meier or the Nelson-Aalen curves, as their
starting point. When these curves overlap early on, only to diverge later,
this will be taken to show delayed effects. When the curves begin to separate
well only to later come back together, this will be taken to indicate a waning
or a reversal of effects. We would not dispute that, but we do feel that the
estimated survival curves are not the best visual tool to examine
nonproportionality. This is very clearly shown in the Long et al. (2017) data.
In light of the observed Kaplan-Meier curves, most of us would take the
proportional hazards assumption to be at least a reasonable approximation to
the relation between the survival experience of the two groups. And yet, the
treatment effect and the bridged process show the actual observations are
poorly described by a proportional hazards assumption.

The consequences of violating the proportional hazards assumption in the Long
et al. (2017) study are not minor. The overall conclusion that dabrafenib and
trametinib are better than placebo is correct, but we only managed to arrive
at this conclusion because the actual effects were so strong and the sample
size was quite large. Had the effect been weaker - yet still very strong under
a more accurate model specification - we may have missed the key scientific
fact that dabrafenib and trametinib are better than placebo. This does not
mean that, despite some shortcomings, we succeeded anyway. A key goal,
detecting a treatment effect, was accomplished. But, the importance of this
effect was largely underestimated. The treatment effect is, in fact, if only
for the first 8.5 months, an order of magnitude greater than what we might
have anticipated from working with the wrongly specified proportional hazards
model. And this is only a part of the story. A nonproportional hazards model,
much closer to the observations as indicated by the satisfactory fit (see
Figure 4 D), suggests that, after week 8.5, there is almost no effect of
treatment. This is easily confirmed by looking at the conditional survival
estimates. Figure 8 shows that for all patients who survive until week 8.5 the
remaining survival is almost independent of group status, the Kaplan-Meier
curves all but overlap.

Figure 8:  The estimated relapse-free survival landmarked at the changepoint
of 8.5 months in the Long et al. (2017) example. Conditioning on living to 8.5
months, there was no further indication of a treatment effect.

Continuing our discussion on the Long et al. (2017) example, we might also
note that from week 8.5 toward the study’s closing date - a considerable
percentage of the total study time - there was no treatment effect. Imagining
that as time goes by and more events are recorded, we will increase the
study’s power, which can correspond to a misassumption in the design. If the
model structure is ignored, our power can decrease with time rather than vice
versa. In practice, the strong early effects, i.e., those observed before week
8.5, are becoming increasingly diluted. A mirror image of this situation can
be seen in immunotherapy studies where an initial period may show either weak
or no effect, followed by a steadily growing effect. Both situations mark
significant departures from a proportional hazards assumption and must be
addressed. How ought we respond to this in the design stage? While we cannot
know the future, little can be predicted with high precision, but we may be
able to use similar or related studies analyzed retrospectively - specifically
making use of the nph2ph transform. We can then consider design strategies
that can help us increase the power, or at least not lose too much power, of
the study that we are planning when nonproportional hazards is a real
possibility.

##  9 Acknowledgements

This work was supported by the NIH Grant P30CA008748.

##  10 Data Availability

The data included in the manuscript are based on digitized Kaplan-Meier
survival curves from published clinical trials. The method detailed by Guyot
et al (2012) was used to digitize each trial presented.

## References

  * André et al. (2020) T. André, K. K. Shiu, T. W. Kim, and et. al.  Pembrolizumab in Microsatellite-Instability-High Advanced Colorectal Cancer.  _N Engl J Med_ , 383(23):2207–2218, 12 2020.
  * Cai and Sun (2003) Z. Cai and Y. Sun.  Local linear estimation for time–dependent coefficients in Cox’s regression models.  _Scandinavian Journal of Statistics_ , 30:93–111, 2003\.
  * Chauvel and O’Quigley (2014) C. Chauvel and J. O’Quigley.  Tests for comparing estimated survival functions.  _Biometrika_ , 101(3):535–552, Sep 2014.
  * Chauvel and O’Quigley (2017) C. Chauvel and J. O’Quigley.  Survival model construction guided by fit and predictive strength.  _Biometrics_ , 73(2):483–494, Jul 2017.
  * Choodari-Oskooei et al. (2012) B. Choodari-Oskooei, P. Royston, and M. K. Parmar.  A simulation study of predictive ability measures in a survival model I: explained variation measures.  _Stat Med_ , 31(23):2627–2643, Oct 2012.
  * Chow (1960) G. Chow.  Tests of Equality Between Sets of Coefficients in Two Linear Regressions.  _Econometrica,_ , 28(3):591–605, Jul 1960.
  * Contal and O’Quigley (1999) C. Contal and J. O’Quigley.  An application of changepoint methods in studying the effect of age on survival in breast cancer.  _Computational Statistics I & Data Analysis_, 30(3):253–270, May 1999.
  * Cox (1972) D. Cox.  Regression models and life tables (with discussion).  _Journal of the Royal Statistical Society, Series B_ , 34:187–220, 1972.
  * Eggermont et al. (2016) A. M. Eggermont, V. Chiarion-Sileni, J. J. Grob, R. Dummer, J. D. Wolchok, H. Schmidt, O. Hamid, C. Robert, P. A. Ascierto, J. M. Richards, C. Lebbé, V. Ferraresi, M. Smylie, J. S. Weber, M. Maio, L. Bastholt, L. Mortier, L. Thomas, S. Tahir, A. Hauschild, J. C. Hassel, F. S. Hodi, C. Taitt, V. de Pril, G. de Schaetzen, S. Suciu, and A. Testori.  Prolonged Survival in Stage III Melanoma with Ipilimumab Adjuvant Therapy.  _N Engl J Med_ , 375(19):1845–1855, 11 2016.
  * Flandre and O’Quigley (2019) P. Flandre and J. O’Quigley.  Comparing Kaplan-Meier Curves with Delayed Treatment Effects: Applications in Immunotherapy Trials.  _J. R. Statist. Soc. C_ , 68(4):915–939, Aug 2019\.
  * Gerds et al. (2008) T. A. Gerds, T. Cai, and M. Schumacher.  The performance of risk prediction models.  _Biom J_ , 50(4):457–479, Aug 2008.
  * Guyot et al. (2012) P. Guyot, A. E. Ades, M. J. Ouwens, and N. J. Welton.  Enhanced secondary analysis of survival data: reconstructing the data from published Kaplan-Meier survival curves.  _BMC Med Res Methodol_ , 12:9, Feb 2012.
  * Gönen and Heller (2005) M. Gönen and G. Heller.  Concordance probability and discriminatory power in proportional hazards regression.  _Biometrika_ , 92:965–970, 2005.
  * Harrell et al. (1996) F. E. Harrell, K. L. Lee, and D. B. Mark.  Tutorial in biostatistics. multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.  _Statistics in medicine_ , 15:361––387, 1996.
  * Hastie and Tibshirani (1990) T. Hastie and R. Tibshirani.  Exploring the nature of covariate effects in the proportional hazards model.  _Biometrics_ , 46:1005–1016, 1990.
  * Jonker et al. (2007) D. J. Jonker, C. J. O’Callaghan, C. S. Karapetis, J. R. Zalcberg, D. Tu, H. J. Au, S. R. Berry, M. Krahn, T. Price, R. J. Simes, N. C. Tebbutt, G. van Hazel, R. Wierzbicki, C. Langer, and M. J. Moore.  Cetuximab for the treatment of colorectal cancer.  _N Engl J Med_ , 357(20):2040–2048, Nov 2007\.
  * Long et al. (2017) G. V. Long, A. Hauschild, M. Santinami, V. Atkinson, M. Mandalà, V. Chiarion-Sileni, J. Larkin, M. Nyakas, C. Dutriaux, A. Haydon, C. Robert, L. Mortier, J. Schachter, D. Schadendorf, T. Lesimple, R. Plummer, R. Ji, P. Zhang, B. Mookerjee, J. Legos, R. Kefford, R. Dummer, and J. M. Kirkwood.  Adjuvant Dabrafenib plus Trametinib in Stage III BRAF-Mutated Melanoma.  _N Engl J Med_ , 377(19):1813–1823, 11 2017.
  * Miller and Siegmund (1982) R. Miller and D. Siegmund.  Maximally selected chi square statistics.  _Biometrics_ , 38(4):1011–1016, 1982.
  * Murphy (1993) S. Murphy.  Testing for a time–dependent coefficient in Cox’s regression model.  _Scandinavian journal of Statistics_ , 20:35–50, 1993.
  * Murphy and Sen (1991) S. A. Murphy and P. K. Sen.  Time-dependent coefficients in a cox-type regression model.  _Stochast. Process. Appl._ , 39:153–180, 1991.
  * O’Quigley (2003) J. O’Quigley.  Khmaladze‐type graphical evaluation of the proportional hazards assumption.  _Biometrika_ , 90(3):269––276, 2003.
  * O’Quigley (2021) J. O’Quigley.  _Survival Analysis: Proportional and Non-Proportional Hazards Regression_.  Springer, 2021.
  * O’Quigley and Flandre (1994) J. O’Quigley and P. Flandre.  Predictive capability of proportional hazards regression.  _Proc Natl Acad Sci U S A_ , 91(6):2310–2314, Mar 1994.
  * O’Quigley and Natarajan (2004) J. O’Quigley and L. Natarajan.  Erosion of regression effect in a survival study.  _Biometrics_ , 60(2):344–51, Jun 2004.
  * O’Quigley and Pessione (1991) J. O’Quigley and F. Pessione.  The problem of a covariate-time qualitative interaction in a survival study.  _Biometrics_ , 47(1):101–115, 1991.
  * Satagopan et al. (2017) J. M. Satagopan, A. Iasonos, and J. G. Kanik.  A reconstructed melanoma data set for evaluating differential treatment benefit according to biomarker subgroups.  _Data Brief_ , 12:667–675, Jun 2017.
  * Scheike and Martinussen (2004) T. Scheike and T. Martinussen.  Maximum likelihood estimation for Cox’s regression model under case–cohort sampling.  _Scandinavian Journal of Statistics_ , 31(2):283–293, 2004.
  * Steyerberg et al. (2013) E. W. Steyerberg, K. G. Moons, D. A. van der Windt, J. A. Hayden, P. Perel, S. Schroter, R. D. Riley, H. Hemingway, D. G. Altman, K. Abrams, D. Altman, A. Briggs, N. Brunner, P. Croft, J. Hayden, A. Hingorani, H. Hemingway, P. Kyzas, N. Malats, K. Moons, G. Peat, P. Perel, R. Riley, I. Roberts, W. Sauerbrei, S. Schroter, E. Steyerberg, A. Timmis, and D. van der Windt.  Prognosis Research Strategy (PROGRESS) 3: prognostic model research.  _PLoS Med_ , 10(2):e1001381, 2013.
  * Uno et al. (2011) H. Uno, T. Cai, M. J. Pencina, R. B. D’Agostino, and L. J. Wei.  On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data.  _Stat Med_ , 30(10):1105–1117, May 2011.
  * Xu and O’Quigley (1999) R. Xu and J. O’Quigley.  A R2 type measure of dependence for proportional hazards models.  _Journal of Nonparametric Statistics_ , 12:83–107, 1999\.