# A Stochastic Compound Failure Model for Testing Resilience of Autonomous
Fixed-Wing Aircraft I: Formulation and Simulation

Thelonious Cooper and Sai Ravela
Earth Signals and Systems Group
Massachusetts Institute of Technology
77 Massachusetts Avenue, Cambridge, MA 02139

###### Abstract

This paper presents a Markov chain model to dynamically emulate the effects of
adverse (failure) flight conditions on fixed-wing, autonomous aircraft system
actuators. It implements a PX4 Autopilot flight stack module that perturbs the
attitude control inputs to the plane’s actuator mixer. We apply this approach
in simulation on a fixed-wing autonomous aircraft to test the controller
response to stochastic compound failures on a range of turning radii.
Statistical measures of the differences between target and simulated flight
paths demonstrate that a well-tuned PID controller remains competitive with
adaptive control in a cascading, compound, transient failure regime.

Adaptive Control, UAV Simulation, Failure Modeling

††preprint: APS/123-QED

##  I Introduction

Despite their inadequacies in adverse conditions, low-cost fixed-wing
Uncrewed/autonomous Aircraft Systems (UAS) are standard in environmental
research. Where satellite imagery is too infrequent, strapping sensors to an
inexpensive UAS and flying them is an attractive means for collecting data. In
these missions, aircraft damage and payload losses are substantial when parts
fail, even when the UAS is repairable. For missions such as flying into a
volcanic plume or shallow cumuli ravela13 ; Ravela2018 , errant ash or soot
jamming a servo could mean data loss and thousands of dollars.

In such flights, numerous things can go wrong. Electrical or communication
failures induced by the environment, low-quality servos and cabling, poor
thermal management, poorly strapped payload, and even mismatched wing surfaces
are all culprits. In many situations, the failures are not immediately
catastrophic but instead appear as incipient transient disturbances that
cascade into compound failures, leading to total failure and permanent loss as
too many things go wrong. Cascading failures often turn initial pilot or
operator confusion or ignorance into panic, contributing to a crash.

While redundant controls and sensors could render the aircraft flyable with
degraded performance when a few parts fail, for autonomous or aided flight,
improving autopilot resilience, handling, and recovering from developing
transient in-flight failures would tremendously advance operating low-cost UAS
in adverse environments. However, this is easier said than done. Hierarchical
PID controllers, at the heart of most stock modern hobby/environmental
research uncrewed/autonomous aircraft autopilots px4 , are generally well-
tuned. Still, conventional wisdom is that finely tuned autopilots are not
robust to “incessant hiccups.” For example, we do not expect stock
hierarchical PID controllers to adapt quickly to changing vehicle dynamics or
failing actuators; some adaptive control appears to be needed.

One approach is carefully modeling aircraft dynamics under various failure
conditions and detecting them in flight to trigger an appropriate control law.
This might involve model predictive control with online model parameter tuning
and additional PID coupling to handle model imperfections. However, such
approaches are likely to be very expensive. Adaptive control must be amenable
to low-cost hardware in the envisioned use cases. In this context, exciting
new developments in Retrospective Cost Adaptive Control (RCAC)
goel2020adaptive ; goel2021a ; goel2020a ; lee2021adaptive  suggest tuning PID
gains in a data-driven manner by optimizing a performance index on a
retrospective performance variable. The RCAC solution recursively updates
gains and is easy to implement on the typically low-cost hardware that stock
autopilots, such as PX4 px4 ; px4reference , used in this research, execute
on.

Figure 1: The Earth Signals and Systems group flies multiple aircraft for
environmental applications, testing new control laws. Here, we show snapshots
of a flight test with RCAC lee2021adaptive . From left to right, the Pixhawk
autopilot mount, the aircraft, ground station, flight path, pitot tube
assembly, and broken nose wheel. Such conditions as taped pitot tubes and
broken assemblies are common in hobby flying.

Previously, we implemented RCAC for quadrotor goel2020adaptive ; goel2021a ;
spencer2021a  and fixed-wing aircraft lee2021adaptive  in the PX4 environment
(PX4 v1.10) px4 ; px4reference  with Holybro Pixhawk (FMUv5) px4HolybroPixhawk
hardware, and have tested them in software-in-the-loop, hardware-in-the-loop,
and assisted flight configurations (see Figure 1). In these tests, we flew the
aircraft in clean and adverse (simulated failure) conditions, where cutting an
actuator simulates failure. These configurations tested single-step failure
injection. Based on our flight testing experience with fixed-wing aircraft, we
posit that in many real-world cases, cascading, compounding intermittent
failures that, left unresolved, become total failures are common.

The current storyline starts with a well-configured aircraft and a well-tuned
autopilot. That is, we are generally satisfied with the autopilot in clean
conditions. We augment this autopilot with an adaptive controller (RCAC) so
that the system shows some skill at adapting and maintaining performance in a
demanding failure regime that is essential to test before use. In this regime,
transient actuator failures emerge, typically milli-seconds to seconds long.

Despite existing solutions for simulating sensor failures, systematic
approaches to testing cascading and compounding actuator failures are not
readily available for PX4-like stock autopilots in small UAS settings. Here,
we develop a stochastic model wherein a stationary Markov chain simulates
actuator failure. It transitions the aircraft from a safe state (the clean
configuration) to various degrees of simulated failure conditions, locking and
releasing a combination of actuators at any given time. In this first attempt,
the Markov transition matrix only loosely corresponds to flight observations
from extensive prior flight experience (for example, see Figure 1). It uses a
Hardware-in-the-Loop (HITL) environment to test simulated failures and
trajectories for fixed-wing aircraft.

Our central hypothesis tests whether RCAC significantly improves resilience
over well-tuned stock autopilots in the transient, cascading, and compound
failure regimes. Our results indicate that with specific hyper-parameter
tuning, RCAC deftly handles step changes to flight configurations. Previous
work lee2021adaptive ; spencer2021a  shows that mixing its control inputs with
the stock PX4 autopilot px4  has a particular advantage for compensating
poorly calibrated PID parameters. However, we find no specific benefit in the
transient failure regime at typical timescales of interest. The adaptations
tightly tune the aircraft to specific regimes, making transitions in a
cascading failure regime difficult. It is possible that the issue could be
addressed by dynamic RCAC hyper-parameter tuning. However, no clear strategies
have emerged as of this writing. Adapting to cascading failures arriving at
different rates for different lengths induces a performance variance that
degrades overall resilience.

The remainder of this paper is as follows. Section II describes related work.
Section III describes the stock PID (PX4 v1.10) px4  controller and RCAC
lee2021adaptive , and Section IV describes flight-path generation, Markov
failure model, and comparison measures. Section V discusses implementation
details for the Holybro Pixhawk (FMUv5), and a Discussion and Conclusions
follow.

##  II Related Work

There several existing adaptive control techniques for controlling fixed-wing
aircraft nguyen2006a  such as sliding mode yu2020a ; fu2021a , backstepping-
based nonlinear controller hirano2019a , and MRAC-based control augmentation
xiong2020a . Although approaches such as ensemble control
trautner2020informative  on embedded GPUs are advancing, computationally
expensive detailed models are still needed. We use RCAC, which is low cost and
uses the past measured data and control inputs to recursively estimate
controller gains rahman2017a ; goel2021a ; lee2021adaptive ; spencer2021a .
Extension to digital PID control kamaldar2019a  and muticopter autopilot
goel2021a ; spencer2021a  are available with an implementation for the PX4
Pixhawk platform. In previous work lee2021adaptive , RCAC was used to
compensate for singleton actuator faults in both calibrated and degraded PX4
conditions lee2021adaptive , which we refer the reader to as the primary
reference for details. Here, we use this framework to test performance over a
stream of the transient compound, cascaded actuator faults.

Prochazka et al. prochazka_development_2021  develop a HITL demonstrator for
validating fault-tolerant control methods for a hybrid UAV, where they
simulate actuator failure and subsequent reconfiguration by a fault-tolerant
flight control law. Jha et al. jha18  develop a fault injector for autonomous
vehicles, including data, hardware, timing, and machine learning subsystem
components. Wen et al. wen19  describe a fault simulator for UAS that includes
sensor and actuator failures with human-in-the-loop flying. Gong et al. gong19
argue that prior data is limited and develop a fault injector for UAS process
health monitoring. In contrast, this work produces a stochastic Markov model
for singleton, compound, and cascading actuator failures. Although
unimplemented here, flight data can train the pre-specified transition
probabilities.

##  III Controllers

Simulation experiments test two flight controllers, the first is the stock PX4
autopilot based on hierarchical PID px4 , and the second is RCAC
lee2021adaptive . This section briefly discusses each controller, and a GitHub
page mitonrgit  contains the implementation of our work in the PX4
environment.

###  III.1 PX4 PID Controllers for Fixed Wing Aircraft

The PX4 autopilot px4reference  receives position setpoints from a Ground
Controller (e.g., QGroundControl, Mission Planner) and controls the fixed-wing
aircraft with two nested-loop cascaded controllers. A position controller
operates the outer loop for longitudinal and lateral motion control, and an
attitude controller operates the inner loop. The position controller’s
longitudinal control section sets the thrust ($T_{s}$) and pitch
($\theta_{s}$) set points, and the lateral control section sets the roll
($\phi_{s}$) set point. The attitude controller cascades two controllers. The
first controller sets pitch ($\dot{\theta}_{s}$) and roll ($\dot{\phi}_{s}$)
rates proportionately to roll and pitch errors (computing the yaw rate is
algebraic). The second controller employs a feedforward and Proportional-
Integral (PI) control law to set the angular acceleration. Control allocation
methods then set actuator deflections using the angular acceleration set
points.

Thus, the position controller takes the True airspeed ($V_{tas,s}$), attitude,
and position setpoints and the corresponding measurements of true airspeed
($V_{tas}$), attitude, and position to set the thrust $T_{s}$, roll
$\phi_{s}$, and pitch $\theta_{s}$ demands. Its longitudinal control section
employs a Total Energy Control System (TECS) bruce1986a ; faleiro1999a ;
lambregts2013a ; argyle2016a , using a different guidance law park2004a  for
lateral control.

The attitude controller’s proportion control law uses gains $k_{\theta}$ and
$k_{\phi}$ to produce rate setpoints:

$\displaystyle\dot{\theta}_{s}$ $\displaystyle=$ $\displaystyle
k_{\theta}(\theta_{s}-\theta_{m}),$ (1) $\displaystyle\dot{\phi}_{s}$
$\displaystyle=$ $\displaystyle k_{\phi}(\phi_{s}-\phi_{m}).$ (2)

The yaw rate setpoint ($\dot{\psi}_{s}$) uses the true airspeed $V_{tas}$ and
attitude to set develop a coordinated turn:

$\dot{\psi}_{s}=\frac{g\tan{\phi_{s}}\cos{\theta_{s}}}{V_{tas}},$ (3)

where $g$ is gravity. The PI control law generates angular acceleration
setpoint $\alpha_{s}$ as:

$\alpha_{s}=\frac{\bar{V}_{tas}}{V_{tas}}\;k_{\omega,\text{ff}}\;\omega_{s}+\left(\frac{\bar{V}_{ias}}{V_{ias}}\right)^{2}G_{\omega,\text{PI}}(\mathbf{q})(\omega_{s}-\omega_{m}),$
(4)

where $\bar{V}_{tas}$ and $\bar{V}_{ias}$ are the trim true and indicated
airspeeds, respectively, $k_{\omega,\text{ff}}$ is the feedforward
proportional gain, and $\omega_{s}$ and $\omega_{m}$ are the desired and
measured body-frame angular velocities, given by

$\omega_{s}=\begin{bmatrix}1&0&sin\theta\\\ 0&cos\phi&sin\phi\;cos\theta\\\
0&-sin\phi&cos\phi\;cos\theta\end{bmatrix}\begin{bmatrix}\dot{\phi}_{s}\\\
\dot{\theta}_{s}\\\ \dot{\psi}_{s}\end{bmatrix}.$ (5)

The PI control law is:

$G_{\omega,\text{PI}}(\mathbf{q})=k_{\omega,P}+{k_{\omega,I}\over\mathbf{q}-1},$
(6)

where $\mathbf{q}$ is the forward-shift operator. The gains $k_{\theta}$,
$k_{\phi}$, $k_{\omega,\text{ff}}$, $k_{\omega,P}$, and $k_{\omega,I}$
constitute 11 variables that must be tuned for flight.

###  III.2 Retrospective Cost Adaptive Control for Fixed-wing Aircraft

Retrospective Cost Adaptive Control lee2021adaptive  is a technique that uses
a retrospective cost function to update the gains of the PID controller
recursively. In the present application, we modify Equations 1 and 2 as:

$\displaystyle\dot{\theta}_{s}$ $\displaystyle=$ $\displaystyle
k_{\theta}(\theta_{s}-\theta_{m})+u_{\theta},$ (7)
$\displaystyle\dot{\phi}_{s}$ $\displaystyle=$ $\displaystyle
k_{\phi}(\phi_{s}-\phi_{m})+u_{\phi}.$ (8)

We modify Equation 4 as :

$\alpha_{s}=\frac{\bar{V}_{tas}}{V_{tas}}\;k_{\omega,\text{ff}}\;\omega_{s}+\left(\frac{\bar{V}_{ias}}{V_{ias}}\right)^{2}G_{\omega,\text{PI}}(\mathbf{q})(\omega_{s}-\omega_{m})+u_{\omega,PI}.$
(9)

The control input terms $u_{\theta}$, $u_{\phi}$ and $u_{\omega,PI}$ are
provided by the RCAC controller.

###  III.3 Adaptive Control Law

To understand how RCAC works, consider the example of a SISO PID controller
with a feedforward term

$\displaystyle u_{k}$ $\displaystyle=$ $\displaystyle
K_{p,k}e_{k-1}+K_{i,k}\nu_{k-1}$ (10)
$\displaystyle+K_{d,k}(e_{k-1}-e_{k-2})+K_{ff,k}r_{k},$

where $k\geq 0$, $e_{k}$ is the error variable, $r_{k}$ is the feed-forward
variable, and

$\displaystyle\nu_{k}$ $\displaystyle\overset{\text{\tiny def}}{=}$
$\displaystyle\sum_{n=0}^{k}g(e_{n})$ (11) $\displaystyle=$
$\displaystyle\nu_{k-1}+g(e_{k-1}).$ (12)

A regression synthesizes the control law:

$u_{k}=\gamma_{k}\rho_{k},$ (13)

where

$\gamma_{k}\overset{\text{\tiny
def}}{=}\begin{bmatrix}e_{k-1}&\nu_{k-1}&(e_{k-1}-e_{k-2})&r_{k}\end{bmatrix},$
(14)

and

$\rho_{k}\overset{\text{\tiny
def}}{=}\begin{bmatrix}K_{p,k}&K_{i,k}&K_{d,k}&K_{ff,k}\end{bmatrix}^{T}.$
(15)

To determine the controller gains $\rho_{k}$, consider a retrospective
performance variable and error model lee2021adaptive

$\hat{e}_{k}(\rho)\overset{\text{\tiny def}}{=}e_{k}+\sigma(\gamma_{k-1}\rho-
u_{k-1}),$ (16)

and the retrospective cost function $J_{k}$

$\displaystyle J_{k}(\rho)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{k}\hat{e}_{n}(\rho)^{T}R_{e}\hat{e}_{n}(\rho)$ (17)
$\displaystyle+(\gamma_{k}\rho)^{T}R_{u}(\gamma_{k}\rho)$
$\displaystyle+(\rho-\rho_{0})^{T}P_{0}^{-1}(\rho-\rho_{0}).$

RCAC seeks the minimizer for $J_{k}$, that is

$\rho_{k+1}\overset{\text{\tiny
def}}{=}\operatorname*{arg\,min}_{\rho}J_{k}(\rho).$ (18)

For $k\geq 0$, $\rho_{k+1}$ is recursively calculated islam2019a

$\displaystyle\rho_{k+1}=\rho_{k}$ (19) $\displaystyle-\sigma
P_{k+1}\gamma_{k-1}^{T}R_{e}[e_{k}+\sigma(\gamma_{k-1}\rho_{k}-u_{k-1})]$
$\displaystyle-P_{k+1}\gamma^{T}_{k}R_{u}\gamma_{k}\rho_{k},$

where

$P_{k+1}=P_{k}-P_{k}\Gamma^{T}_{k}(R_{ue}^{-1}+\Gamma_{k}P_{k}\Gamma_{k}^{T})^{-1}\Gamma_{k},$
(20)

and

$\Gamma_{k}\overset{\text{\tiny
def}}{=}\begin{bmatrix}\sigma\gamma_{k-1}&\rho_{k}\end{bmatrix}^{T},$ (21)

and

$R_{ue}\overset{\text{\tiny def}}{=}\begin{bmatrix}R_{e}&0\\\
0&R_{u}\end{bmatrix}.$ (22)

The following equation provides the control input

$u_{k+1}=\gamma_{k+1}\rho_{k+1}$ (23)

The SISO case extends to MIMO through various parameterizations goel2020a .
Here, RCAC modifies the PX4 attitude controller as shown in Equations 7 \- 9,
where the control inputs $u_{\theta}$, $u_{\phi}$, and $u_{\omega,PI}$ are
generated for every discrete time step. We refer the reader lee2021adaptive
for additional details of RCAC. The following sections discuss details of
applying the stock PX4 and RCAC to flight tests.

##  IV Methods

Testing the stock and adaptive regimes requires three additional components.
They include space-filling curves at different scales to express flight plans,
the stochastic compound failure model, and a measure for comparing the desired
and executed flight paths. These are discussed in this section.

###  IV.1 Hilbert Curves as Flight Test Patterns

A Hilbert curve hilbert1891  is a continuous fractal space-filling curve,
typically constructed as a limit of piecewise lines. The Hausdorff dimension
of Hilbert curves is two, and its graph is a compact set homeomorphic to the
closed unit interval, with Hausdorff dimension 2. The length of the $n^{th}$
curve is $2^{n}-{1\over 2^{n}}$, i.e., the size grows exponentially with $n$,
but each curve is in the unit square wikipediaHilbertCurve . The following
pseudo code mathworksHilbertCurve  is an easy example of Hilbert curves:

$[x,y]\leftarrow Hilbert(n)$

if $n\leq 0$ then

$x\leftarrow 0$

$y\leftarrow 0$

else

$[xx,yy]\leftarrow Hilbert(n-1)$

$x\leftarrow\frac{1}{2}*[-\frac{1}{2}+yy\;\;-\frac{1}{2}+xx\;\;\frac{1}{2}+xx\;\;\frac{1}{2}-yy]$

$y\leftarrow\frac{1}{2}*[-\frac{1}{2}+xx\;\;\frac{1}{2}+yy\;\;\frac{1}{2}+yy\;\;-\frac{1}{2}-xx]$

end if

We designed a flight plan to test various turning frequencies based on Hilbert
paths of varying order. We implemented the generation of these curves and used
curves of orders 1,2,3 and 4 for each quadrant, respectively, see Figures 6
and 7.

Figure 2: Failure states and the corresponding actuator failures. Each row
depicts a Boolean condition indicating simulated failure (1) “Gum” of
particular actuators or combinations. An $X$ indicates that the row (failure)
is never simulated (NA). The failure simulator proceeds if the handheld switch
(SW on channel 8 in our setup) is on.  Figure 3: Markov process transition
diagram, downward arrows represent edges to the ground state. The state
numbers correspond to the table in Figure 2.

###  IV.2 Simulating Compound Failures

The compound failure model tests multiple combinations of failures using a
Markov chain. The ground state of No Failure: 0 is continued indefinitely
until a switch (S/W) on the ground handheld controller is activated to
initiate failures. Turning the switch off reverts the aircraft to $State-0$
while continuing the flight.

Each transition from the ground state induces a single actuator failure.
Although many forms of failure are admissible through our design, we only
implemented the corresponding actuator getting stuck (“gummed”) aileron left
(AIL-L), right (AIL-R), an elevator (ELE), throttle (THR) or Rudder (Rud).
Thus, the corresponding actuator is gummed up once transitioned to these
states from the ground state. The aircraft may stay in this state, revert to
State 0, or transition to compound failure states with one more failure. After
that, transitions can continue through various combinations across states,
each additional transition adding or removing a failure, in addition to
staying in the same state or returning to the ground state. Certain compound
states are never allowed. For example, State 10 in Figure 2, where both the
elevator and throttle are stuck, will render the aircraft unrecoverable
(depending on the throttle setting). Similarly, for field safety reasons, both
ailerons are not failed (though one could still fly with the rudder in some
cases). Thus, all the states are in some sense “flyable” while still testing
failures.

Our Markov chain simulations used a loop-back probability of $0.3$ and a
ground-state transition probability of $0.4$, distributing the transition
edges to valid additional failure states equally amongst the remaining
probability. The transitions are discrete in time and triggered every second.
Thus, the aircraft must stay in a state for at least one second with a $30\%$
chance of remaining there in the next second, $60\%$ possibility of
transitioning to another failure, or a $40\%$ chance of returning to the
ground state. In this way, not only is there persistence of a single failure,
but various compound and intermittent failures are also triggered.

Ensemble simulations of stochastic adverse flight conditions yield performance
statistics richer than single-perturbation tests (e.g., lee2021adaptive ) for
characterizing performance. In comparing different control schemes, each
method can be subject to the same failure sequence to assess sample-wise and
ensemble-wise differences.

###  IV.3 Comparing Flight Paths Using Dynamic Time Warping

To compare flight paths, we resample the reference flight plan to the same
length as the flown trajectory, saving it at nearly constant time intervals.
Since our objective is to assess an ensemble of flight trajectories, the
existence of relative timing differences can lead to the overall inflation of
covariance ravela2015stics . The Kantorovich-Rubenstein or Wasserstein metric
kantorovich , offers a natural approach to dealing with timing errors arriving
at features (e.g., sharp turns). Here, we analyze the difference between the
aircraft’s performance in regular flight vs. facing a stochastic failure
process using Dynamic Time Warping(DTW) dtw01  between the flight paths. DTW
was first proposed within the controls community in the 60s bellman59  and has
since enjoyed broad popularity in comparing discrete signals with index
alignment errors. Note that DTW is not a metric, and in that, it differs from
the Wasserstein distance. However, it is quick to compute, the values are
usually close, and the reference signal (flight plan) is the same. A low DTW
distance between the affected and unaffected flight paths represents that the
autopilot could maintain course well despite the actuator failures. A high DTW
distance implies the opposite. The nominal algorithm for DTW is as follows
wikipediaDynamicTime :

DTWfunc( $x[1\ldots m]$, $y[1\ldots n]$)

$DTW[i,j]=\infty,\;\;\forall i=1\ldots m,j=1\ldots n$

$DTW[0,0]=0;$

for i=1:m do

for j=1:n do

$DTW[i,j]=Distance(x[i],y[j])+\min(DTW[i-1,j],DTW[i,j-1],DTW[i-1,j-1])$

end for

end for

The function $Distance$ is a nominal point-wise Euclidean distance calculation
or a metric in Hilbert space. The all-pair calculations incur quadratic costs,
and incorporating window/taper functions simplifies it.

##  V Implementation

Figure 4: Modified PX4 architecture showing the failure injection module (in
orange).

To implement the control perturbation, we created a new instance of the
actuator_controls channel in the publish-subscribe px4 model. We then added a
vehicle command that would tell the simulator or mixer to listen to the new
instance of the actuator_controls topic. Our module then listens to the
original actuator_controls topic and republishes it or sends a zero indicating
a broken link. Besides creating our module, we had to modify the px4
architecture by altering the simulator module to listen for a command to
switch control topic instances.

Figure 5: The PX4 implementation px4reference  on PixHawk px4HolybroPixhawk
splits the computation between two processors. Failure streams were re-
implemented on FMU because the IO processor was weaker and could not sustain
the computation.

The Pixhawk PX4 (model) we tested with hardware in the loop testing uses two
processors to perform flight control. The primary processor handles all of the
mathematical computation, and a weaker coprocessor handles output operations
and relays sensor information to the primary processor. We initially tried
implementing the failure simulation on the IO coprocessor. However, this
attempt failed because the additional processing of the failure simulation
caused the coprocessor to break its real-time condition and crash. The IO
coprocessor operates by a serial interface that blind writes into hard-coded
memory sections which a critical section reads before being sent to actuators
and other output devices.

One challenge we encountered and could not solve fully was a race condition
within the DDS data bus when stopping the failure simulation. Due to the real-
time constraints of the system, after instructing the failure simulation
process to stop, clock cycles are offered for a limited amount of time before
deallocating its memory. This time is, unfortunately, faster than the update
rate of the simulator, which waits for a signal to switch back to the original
control signal instance only after attempting to read from the deallocated
phony control signal buffer. That meant we could not stop the failure
simulation process once it started. We could, however, command it not to
modify the control inputs it parses and pass them along to the mixer
unmodified, effectively disabling the purpose of the process.

##  VI Experiments

Figure 6: Examples of Nominal Flight Paths for PX4 (left) and RCAC (right).
Figure 7: Examples of Failure Simulated Flight Paths for PX4 (left) and RCAC
(right).

The PX4 autopilot is well-tuned for many hobby UAS. We calculated these gains
for the aircraft under consideration 1 through extensive flight testing
starting from readily available gains for similar aircraft. For RCAC, we
performed several flight tests. We checked whether the parameter led to
unstable flight in assisted mode with random flight-stick inputs for each
controlled attitude variable. In many cases, the onset of instability occurred
a considerable time after the stick perturbation. Therefore, we set nominal
values, which may slow the gain adaptation and be sub-optimal. Table 3 lists
the settings we used.

PTCH | P0 | 1
---|---|---
PTCH | RU | 0.001
PITCH RATE | P0 | 1e-4
PITCH RATE | RU | 0.1
ROLL | P0 | 1
ROLL | RU | 0.001
ROL RATE | P0 | 1e-4
ROLL RATE | RU | 0.1
YAW RATE | P0 | 1e-4
YAW RATE | RU | 0.1
Table 1: Nominal RCAC (hyper) parameters. The mixing variable $\alpha$ is to
$1$ (or $0.5$).

For the standard PX4 PID-based autopilot, we ran $8$ trials of the stochastic
adversarial actuator failure process. For an experimental Retrospective Cost
Adaptive Control(RCAC) algorithm, we ran $11$ trials with the mixing
hyperparameter $\alpha=0.5$ and 12 trials with $\alpha=1$. As a baseline, we
also ran several clean flights of the PX4 PID controller. We saw they had an
average self-similarity under the DTW metric of $133$ meters with a standard
deviation of $47$ meters. We then did the same for the RCAC $\alpha=0.5$
autopilot and saw an average self-similarity under the DTW metric of $172$
meters with a standard deviation of $46$ meters. And then again did the same
for the RCAC $\alpha=1$ autopilot and saw that they had an average self-
similarity under the DTW metric of $136$ meters with a standard deviation of
$35$ meters.

The results of our experiments are summarized in the following tables

Flight Configuration | Mean | Standard Deviation
---|---|---
Px4 PID | 461 | 39
RCAC $\alpha=0.5$ | 604 | 45
RCAC $\alpha=1$ | 457 | 19
Table 2: Mean and standard dev of DTW distance to target flight path for unperturbed autopilot configurations Flight Configuration | Mean | Standard Deviation
---|---|---
Px4 PID | 595 | 154
RCAC $\alpha=0.5$ | 1196 | 302
RCAC $\alpha=1$ | 610 | 163
Table 3: Mean and standard dev of DTW distance to target flight path for
failure-simulated autopilot configurations Figure 8: Box plot of DTW distances
for nominal (clean) and failure simulated (adverse) flight conditions.

##  VII Conclusions and Future Work

This paper considers transient singleton, compound, and cascading actuator
failure scenarios and tests the performance of a stock autopilot PX4 (v1.10)
simulation for Holybro Pixhawk (FMUv5) hardware and fixed-wing aircraft. In
actual flight experiments we have conducted elsewhere, multiple failures often
cascade. For example, transient spikes in electrical load can trigger losses
anywhere in the system. The PX4 autopilot, like others, is often tuned for the
aircraft, and considerable attention is devoted to this effort before
commencing flight missions. The case for extending RCAC-like adaptive control
from tuning degraded gains to failure streams emerging in well-tuned
autopilots is interesting and essential. The operational setting is of a tuned
autopilot with the RCAC as “safety” that would compensate for the transient
failure stream.

In this paper, a Markov chain simulates failures with a stationary
distribution that reaches all clean and adverse modes. The Markov chain is
designed not to cause total catastrophic failure of the aircraft. Simulated
Hilbert curves are flight paths for testing PX4 and PX4-RCAC, with mixing
parameter $\alpha=0.5$ and $\alpha=1$. The PX4 flight parameters were
determined for the model aircraft through flight testing and set in
simulations. We selected RCAC hyerparameters that led to stable flight over
long duration for the same model aircraft. DTW was used to compare resampled
flight plans and simulated flight paths to account for timing-error artifacts
between ensemble simulations. The ensemble simulations calculate the mean
position error and the error variance. Under clean configurations, RCAC
slightly outperforms the stock PX4 PID controller, but under the hostile
regime, it under performs. Additional, possibly dynamic, hyper-parameter
tuning may be required to deal with intermittent failure streams. With the
extra work, we will implement this in actual flight tests in the future. We
also anticipate learning the Markov chain using flight logs gleaned from the
PX4 community.

##  VIII Contributions and Acknowledgment

Corresponding Author: Thelonious Cooper (theloni@mit.edu). Thelonious Cooper
wrote code and performed the experiments, Sai Ravela developed the transition
tables, and Erina Yamaguchi reported the original Hilbert path code. Cooper
and Ravela drafted and edited the paper. We thank Prof. Ankit Goel, Jyonghoon
Lee, Juan Augusto Paredes, John Spencer, and Prof. Dennis Bernstein for their
support and collaboration. We are grateful for funding from ONR
(N00014-19-1-2273) and the MIT Weather Extreme and CREWSNET Climate Grand
Challenge projects.

## References

  * (1) Dynamic time warping - Wikipedia — en.wikipedia.org.  https://en.wikipedia.org/wiki/Dynamic_time_warping.  [Accessed 14-May-2023].
  * (2) Hilbert Curve — mathworks.com.  https://www.mathworks.com/matlabcentral/fileexchange/4646-hilbert-curve.  [Accessed 14-May-2023].
  * (3) Hilbert curve - Wikipedia — en.wikipedia.org.  https://en.wikipedia.org/wiki/Hilbert_curve.  [Accessed 14-May-2023].
  * (4) Holybro Pixhawk 4 (FMUv5): PX4 User Guide.  https://docs.px4.io/v1.11/en/flight_controller/pixhawk4.html.  [Accessed 14-May-2023].
  * (5) Dynamic Time Warping, pages 69–84.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2007.
  * (6) PX4 Reference Flight Controller Design, PX4 User Guide.  https://docs.px4.io/main/en/hardware/reference_design.html, 2022\.  [Accessed 14-May-2023].
  * (7) M.E. Argyle and R.W. Beard.  Nonlinear total energy control for the longitudinal dynamics of an aircraft.  In 2016 American Control Conference (ACC, page 6741–6746.
  * (8) R. Bellman and R. Kalaba.  On adaptive control processes.  IRE Transactions on Automatic Control, 4(2):1–9, 1959.
  * (9) K. Bruce, J. Kelly, and J. Person.  Nasa b737 flight test results of the total energy control system.  In Astrodynamics Conference, page 2143.
  * (10) Thelonious Cooper, Erina Yamaguchi, and Sai Ravela.  ESSG (MIT) PX4 Failure Simulator Flight Stack.  https://github.com/theloni-monk/PX4-Autopilot-ESSG-FWFA, 2022.  [Accessed 14-May-2023].
  * (11) L. Faleiro and A. Lambregts.  Analysis and tuning of a total energy control system control law using eigenstructure assignment.  Aerospace science and technology, 3(3):127–140,.
  * (12) M. Fu, Z. Yu, and Y. Zhang.  Adaptive fault-tolerant control of fixedwing uav under actuator saturation and state constraints.  In 2021 International Conference on Unmanned Aircraft Systems (ICUAS, page 47–52.
  * (13) A. Goel, S.A.U. Islam, and D.S. Bernstein.  Adaptive control of mimo systems using sparsely parameterized controllers.  In 2020 American Control Conference (ACC), page 5340–5345.
  * (14) A. Goel, J.A. Paredes, H. Dadhaniya, S.A.Ul Islam, A.M. Salim, S. Ravela, and D. Bernstein.  Experimental implementation of an adaptive digital autopilot.  In 2021 American Control Conference (ACC), page 3737–3742.
  * (15) Ankit Goel, Abdulazeez Mohammed Salim, Ahmad Ansari, Sai Ravela, and Dennis Bernstein.  Adaptive digital pid control of a quadcopter with unknown dynamics.  arXiv:2006.00416, 2020.
  * (16) Siyang Gong, Shengwei Meng, Benkuan Wang, and Datong Liu.  Hardware-in-the-loop simulation of uav for fault injection.  In 2019 Prognostics and System Health Management Conference (PHM-Qingdao), pages 1–6, 2019.
  * (17) D. Hilbert.  Über die stetige abbildung einer linie auf ein flächenstück.  Mathematische Annalen, (38):459–460, 1891.
  * (18) S. Hirano, K. Uchiyama, and K. Masuda.  Controller design using backstepping algorithm for fixed-wing uav with thrust vectoring system.  In 2019 International Conference on Unmanned Aircraft Systems (ICUAS, page 1084–1088. IEEE.
  * (19) S.A.U. Islam and D.S. Bernstein.  Recursive least squares for real-time implementation.  IEEE Control Systems Magazine, 39(3):82–85,.
  * (20) Saurabh Jha, Subho S. Banerjee, James Cyriac, Zbigniew T. Kalbarczyk, and Ravishankar K. Iyer.  Avfi: Fault injection for autonomous vehicles.  In 2018 48th Annual IEEE/IFIP International Conference on Dependable Systems and Networks Workshops (DSN-W), pages 55–56, 2018.
  * (21) M. Kamaldar, S.A.U. Islam, S. Sanjeevini, A. Goel, J.B. Hoagg, and D.S. Bernstein.  Adaptive digital pid control of first-order-lag plus-dead-time dynamics with sensor, actuator, and feedback nonlinearities.  Advanced Control for Applications, 1(1):20,.
  * (22) L. V. Kantorovich.  Mathematical methods of organizing and planning production.  Management Science, 6(4):366–422, 1960.
  * (23) A.A. Lambregts.  Tecs generalized airplane control system design–an update.  In Advances in Aerospace Guidance, Navigation and Control, page 503–534. Springer.
  * (24) Joonghyun Lee, John Spencer, Juan Augusto Paredes, Sai Ravela, Dennis S. Bernstein, and Ankit Goel.  An adaptive digital autopilot for fixed-wing aircraft with actuator faults.  arXiv:2110.11390, 2021.
  * (25) Lorenz Meier, Dominik Honegger, and Marc Pollefeys.  Px4: A node-based multithreaded open source robotics framework for deeply embedded platforms.  In 2015 IEEE International Conference on Robotics and Automation (ICRA), pages 6235–6240, 2015.
  * (26) N. Nguyen, K. Krishnakumar, J. Kaneshige, and P. Nespeca.  Dynamics and adaptive control for stability recovery of damaged asymmetric aircraft.  In AIAA Guidance, navigation, and control Conference and Exhibit, page 6049.
  * (27) S. Park, J. Deyst, and J. How.  A new nonlinear guidance logic for trajectory tracking.  In AIAA guidance, navigation, and control conference and exhibit, page 4900.
  * (28) F. Prochazka, S. Krüger, G. Stomberg, and M. Bauer.  Development of a hardware-in-the-loop demonstrator for the validation of fault-tolerant control methods for a hybrid UAV.  CEAS Aeronautical Journal, 12(3):549–558, August 2021.
  * (29) Y. Rahman, A. Xie, and D.S. Bernstein.  Retrospective cost adaptive control: Pole placement, frequency response, and connections with lqg control.  IEEE Control System Magazine, 37(5):28–69,.
  * (30) Sai Ravela.  Mapping coherent atmospheric structures with small unmanned aircraft systems.  In AIAA Infotech@Aerospace (I@A) Conference, pages 1–11, 2013.
  * (31) Sai Ravela.  Statistical Inference for Coherent Fluids.  LNCS, 8964:121–133, 2015.
  * (32) Sai Ravela.  Tractable Non-Gaussian Representations in Dynamic Data Driven Coherent Fluid Mapping, pages 29–46.  Springer International Publishing, Cham, 2018.
  * (33) J. Spencer, J. Lee, J.A. Paredes, A. Goel, and D. Bernstein.  An adaptive pid autotuner for multicopters with experimental results.  arXiv:2109.12797.
  * (34) Margaret Trautner, Gabriel Margolis, and Sai Ravela.  Informative neural ensemble kalman learning.  arXiv:2008.09915, 2020.
  * (35) Jiayun Wen, Hongxia Ji, Honglun Wang, Menghua Zhang, Dawei Li, and Jianfa Wu.  Design of a real-time uav fault injection simulation system.  In 2019 IEEE International Conference on Unmanned Systems (ICUS), pages 767–772, 2019.
  * (36) J. Xiong, Y. Yang, Z. Cheng, L. Liu, Y. Wang, and H. Fan.  Observer like model reference adaptive augmenting based fixed-wing uav control.  In 2020 39th Chinese Control Conference (CCC), page 6804–6809.
  * (37) Z. Yu, H. Badihi, Y. Zhang, Y. Ma, B. Jiang, and C.-Y. Su.  Fractional order sliding-mode fault-tolerant neural adaptive control of fixed-wing uav with prescribed tracking performance.  In 2020 2nd International Conference on Industrial Artificial Intelligence (IAI), page 1–6.


# Bayesian Spatial Homogeneity Pursuit for Survival Data with an Application
to the SEER Respiratory Cancer Data

Lijiang Geng
Department of Statistics, University of Connecticut, CT 06269, USA  Guanyu Hu
Department of Statistics, University of Missouri - Columbia, MO 65211, USA

###### Abstract

In this work, we propose a new Bayesian spatial homogeneity pursuit method for
survival data under the proportional hazards model to detect spatially
clustered patterns in baseline hazard and regression coefficients. Specially,
regression coefficients and baseline hazard are assumed to have spatial
homogeneity pattern over space. To capture such homogeneity, we develop a
geographically weighted Chinese restaurant process prior to simultaneously
estimate coefficients and baseline hazards and their uncertainty measures. An
efficient Markov chain Monte Carlo (MCMC) algorithm is designed for our
proposed methods. Performance is evaluated using simulated data, and further
applied to a real data analysis of respiratory cancer in the state of
Louisiana.

Keywords: Geographically Weighted Chinese Restaurant Process, MCMC, Piecewise
Constant Baseline Hazard, Spatial Clustering

##  1 Introduction

Clinical data on individuals are often collected from different geographical
regions and then aggregated and analyzed in public health studies. The most
popular dataset is the Surveillance, Epidemiology, and End Results (SEER)
program (SEER, 2016) data which routinely collects population-based cancer
patient data from 20 registries across the United States. This data provides
prognostic and demographic factors of cancer patients. In this paper, we focus
our study on Louisiana respiratory cancer data which was analyzed in Mu et al.
(2020). Analysis of such data conducted on a higher level often assume that
covariate effects are constant over the entire spatial domain. This is a
rather strong assumption, as all intrinsic heterogeneities in data are
ignored. For example, if one was to study the hazard for patients with lung
cancer, it is expected that the true hazard is not the same in areas where
there is little air pollution and severe air pollution, even for patients with
similar characteristics. From Tobler’s first law of geography (Tobler, 1970),
it is reasonable to consider similarities between nearby locations in survival
data due to environmental circumstances in geographically close regions. In
this paper, we will recover the spatial homogeneity pattern of respiratory
cancer survival rates among different counties in state of Louisiana.

Existing approaches that account for such patterns in survival data can be put
into two major categories. The first one is to incorporate spatial random
effects in survival models such as the accelerated failure time (AFT) model
and the proportional hazards model (Banerjee et al., 2003; Banerjee and Dey,
2005; Zhou et al., 2008; Zhang and Lawson, 2011; Henderson et al., 2012), such
that spatial variations are accounted for by different intercepts for
different regions, while parameters for covariates are held constant. Another
important approach, instead of assuming all covariate effects are constant,
allows parameters to be spatially varying in parametric, nonparametric, and
semiparametric models (Hu et al., 2020; Hu and Huffer, 2020; Xue et al.,
2020).

Despite their flexibility, the aforementioned spatially varying coefficients
models can be unnecessarily large. Imposing certain constraints on nearby
regions so that they have the same parameter values provides an efficient way
of reducing the model size without sacrificing too much of its flexibility.
While similar endeavor have been made to cluster spatial survival responses
(Huang et al., 2007; Bhatt and Tiwari, 2014), the clustering of covariate
effects and baseline hazards have yet to be studied for survival data.

Two challenges are to be tackled for clustering of coefficients and baseline
hazards for spatial survival models. First, the spatial structure needs to be
appropriately incorporated into the clustering process. Contiguousness
constraints should be added so that truly similar neighbors are driven to the
same cluster. The constraints, however, should not be overly emphasized, as
two distant regions may still share similar geographical and demographical
characteristics and thus parameters. Existing methods, such as in Lee et al.
(2017, 2019) and Li and Sang (2019), do not allow for globally discontiguous
clusters, which is a serious limitation. Second, the true number of clusters
is unknown, and needs to be estimated. With the probabilistic Bayesian
framework, simultaneous estimation of the number of clusters and the
clustering configuration for each region is achieved by complicated search
algorithms (e.g., reversible jump MCMC, Green, 1995) in variable dimensional
parameter spaces. Such algorithms assign a prior to the number of clusters
that needs to be updated in every MCMC iteration, which made them difficult to
implement or automate, and suffer from mixing issues as well as lack of
scalability. Nonparametric Bayesian approaches, such as the Chinese restaurant
process (CRP; Pitman, 1995), provide another approach to allow for
uncertainties in the number of clusters. Its extension, the distance dependent
CRP (ddCRP; Blei and Frazier, 2011), considers spatial information, and makes
a flexible class of distributions over partitions that allows for dependencies
between their elements. The CRP framework, however, has been shown to be
inconsistent in its estimation of number of clusters (Miller and Harrison,
2013). Lu et al. (2018) proposed the powered CRP that suppresses the small
tail clusters. Similar to the traditional CRP, however, it does not consider
distance information, and therefore is not well-suited when spatial
homogeneity is to be detected.

To address these challenges, in this work, we consider a spatial proportional
hazards model, and propose a geographically weighted Chinese restaurant
process (gwCRP) to capture the spatial homogeneity of both the regression
coefficients and baseline hazards over subareas under piecewise constant
hazards models framework (Friedman et al., 1982). Our main contributions in
this paper are three folds. First, we develop a new nonparametric Bayesian
method for spatial clustering which combines the ideas of geographical weights
and Dirichlet mixture models to leverage geographical information. Compared
with existing methods, our proposed approach is able to capture both locally
spatially contiguous clusters and globally discontiguous clusters. Second, an
efficient Markov chain Monte Carlo (MCMC) algorithm is proposed for our
proposed model without reversible jumps to simultaneously estimate the number
of clusters and clustering configuration. In addition, we apply our method to
the analysis of the Surveillance, Epidemiology, and End Results (SEER) Program
data in the state of Louisiana among different counties, which provide
important information to study spatial survival rates.

The remainder of the paper is organized as follows. In Section 2, we develop a
homogeneity pursuit of survival data in the piecewise constant proportional
hazard framework with gwCRP prior. In Section 3, a collapsed Gibbs sampler
algorithm and post MCMC inference are discussed. The extensive simulation
studies are carried out in Section 4. For illustration, our proposed
methodology is applied to respiratory cancer survival data in Section 5.
Finally, we conclude this paper with a brief discussion in Section 6.

##  2 Methodology

###  2.1 Spatial Piecewise Constant Hazards Models

Let $T_{\ell i}$ denote the survival time for patient $\ell$ at location
$s_{i}$, with $\delta_{\ell i}=1$ representing the event and $\delta_{\ell
i}=0$ indicating censored, and $X_{\ell}(s_{i})$ denotes the vector of
covariates corresponding to $T_{\ell i}$ for $i=1,2,...,n$, and
$\ell=1,2,...,n_{i}$, where $n_{i}$ denotes the number of the patients at
location $s_{i}$. In this paper, $s_{1},s_{2},\ldots,s_{n}$ are areal units
which is defined in Banerjee et al. (2014). Let $\bm{D}=\\{(T_{\ell
i},\delta_{\ell i},X_{\ell}(s_{i})),i=1,2,...,n,\ell=1,2,...,n_{i}\\}$. denote
the observed data. We consider a proportional hazards model (Cox, 1972) with
piecewise constant baseline hazard. We partition $[0,\infty)$ into $J$
intervals ($0=a_{0}<a_{1}<\dots<a_{J}=\infty$), then the hazard function is
given by

$\lambda(t|X_{\ell}(s_{i}))=\lambda_{0}(t)\exp(X_{\ell}(s_{i})^{\top}{\mbox{\boldmath$\beta$}}),\quad$
(1)

with piecewise constant baseline hazard function $\lambda_{0}(t)=\lambda_{j}$
for $a_{j-1}\leq t<a_{j},~{}j=1,\ldots,J$.

For the piecewise constant hazard function mentioned in (1), the baseline
hazards $\lambda_{1},\dots,\lambda_{J}$ and regression coefficients $\beta$
are constants over different regions. Due to observed environmental factors,
spatially varying patterns in baseline hazards and regression coefficients of
hazard function need to be considered. The piecewise constant hazard function
with spatially varying pattern is therefore given by

$\lambda(t|X_{\ell}(s_{i}))=\lambda_{0(s_{i})}(t)\exp(X_{\ell}(s_{i})^{\top}{\mbox{\boldmath$\beta$}}(s_{i})),\quad$
(2)

where $\lambda_{0(s_{i})}(t)=\lambda_{j}(s_{i})$ for $a_{j-1}\leq
t<a_{j},~{}j=1,\ldots,J$. Under this model,
${\mbox{\boldmath$\lambda$}}(s_{i})=(\lambda_{1}(s_{i}),\dots,\lambda_{J}(s_{i}))^{\top}$
and ${\mbox{\boldmath$\beta$}}(s_{i})$ represent the location-specific
baseline hazards and regression coefficients.

After some algebra, the logarithm of likelihood function for observed survival
data $\bm{D}$ is obtained as

$\begin{split}&\log\mathcal{L}({\mbox{\boldmath$\beta$}}(s_{i}),{\mbox{\boldmath$\lambda$}}(s_{i}),i=1,\ldots,n\mid\bm{D})\\\
&=\sum_{i=1}^{n}\left\\{\sum_{j=1}^{J}d_{ji}\log\lambda_{j}(s_{i})+\sum_{\ell=1}^{n_{i}}\delta_{\ell
i}X_{\ell}(s_{i})^{\top}{\mbox{\boldmath$\beta$}}(s_{i})-\sum_{j=1}^{J}\lambda_{j}(s_{i})\left[\sum_{\ell=1}^{n_{i}}\Delta_{j}(T_{\ell
i})\exp(X_{\ell}(s_{i})^{\top}{\mbox{\boldmath$\beta$}}(s_{i}))\right]\right\\},\end{split}$
(3)

where $d_{ji}=\sum_{\ell=1}^{n_{i}}\delta_{\ell
i}{\bf{{1}}}_{[a_{j-1},a_{j})}(T_{\ell i})$, which represents the number of
people at location $s_{i}$ who experience the event during the time period
from $a_{j-1}$ to $a_{j}$, and $\Delta_{j}(t)=t-a_{j-1}\text{ for }a_{j-1}\leq
t<a_{j}.$

For one particular location $s_{i}$, let
$\bm{\eta}(s_{i})=\log{\mbox{\boldmath$\lambda$}}(s_{i})$ and define
${\mbox{\boldmath$\theta$}}(s_{i})=({\mbox{\boldmath$\beta$}}(s_{i})^{\top},\bm{\eta}(s_{i})^{\top})^{\top}$
the collection of parameters, then the maximized likelihood estimate (MLE)
$\widehat{{\mbox{\boldmath$\theta$}}}(s_{i})$ can be obtained by solving the
score function, which is the derivative of the logarithm of likelihood
function in (3), and the estimated variance-covariance matrix of MLE is
$\widehat{\Sigma}_{i}=(-H)^{-1}$, respectively, where $(-H)$ denotes the
negative Hessian matrix. Based on the MLEs and estimated variance-covariance
matrices, we have the following approximation of the likelihood.

######  Proposition 1.

We assume the regularity conditions A-D in Friedman et al. (1982). As
$n_{i}\rightarrow\infty,\,i=1\ldots,n$, the data likelihood
$\mathcal{L}({\mbox{\boldmath$\beta$}}(s_{i}),{\mbox{\boldmath$\lambda$}}(s_{i}),i=1,\ldots,n\mid\bm{D})$
is approximated as

$\mathcal{L}({\mbox{\boldmath$\beta$}}(s_{i}),{\mbox{\boldmath$\lambda$}}(s_{i}),i=1,\ldots,n\mid\bm{D})\approx\prod_{i=1}^{n}\text{MVN}(\widehat{{\mbox{\boldmath$\theta$}}}(s_{i})|{\mbox{\boldmath$\theta$}}(s_{i}),\widehat{\Sigma}_{i}),$
(4)

where MVN stands for the multivariate normal distribution.

The derivations of $\widehat{\Sigma}_{i}$ and the proof of Proposition 1 are
given in Section A and Section B of Supporting information. Instead of using
the log likelihood in (3), our following model is based on normal
approximation in Proposition 1 for computational convenience.

Based on the normal approximation given in Proposition 1, a natural way which
follows Gelfand et al. (2003) for spatially varying pattern of baseline
hazards and regression is to give a Gaussian process prior to
${{\mbox{\boldmath$\theta$}}}(s_{i}),\,i=1,\ldots,n$. The Gaussian process for
${{\mbox{\boldmath$\theta$}}}(s_{i}),\,i=1,\ldots,n$ is defined as

${\mbox{\boldmath$\theta$}}\sim\text{MVN}(\bm{1}_{n\times
1}\otimes\bm{\mu},\bm{H}(\phi)\otimes\Sigma),$ (5)

where
${\mbox{\boldmath$\theta$}}=({\mbox{\boldmath$\theta$}}(s_{1})^{\top},\ldots,{\mbox{\boldmath$\theta$}}(s_{n})^{\top})^{\top}$,
$\bm{\mu}$ is a $p+J$ dimensional vector, $\bm{H}(\phi)$ is a $n\times n$
spatial correlation matrix depending on the distance matrix with parameter
$\phi$, $\Sigma$ is a $(p+J)\times(p+J)$ covariance matrix, and $\otimes$
denotes Kronecker product. The $(i,j)$-th entry of $\bm{H}(\phi)$ is
$\exp(-\phi|s_{i}-s_{j}|)$, where $|s_{i}-s_{j}|$ is the distance between
$s_{i}$ and $s_{j}$, and $\phi>0$ is the range parameter for spatial
correlation. For the Gaussian process prior, the parameters of closer
locations have stronger correlations.

For many spatial survival data, some regions will share same covariate effects
or baseline hazards with their nearby regions. In addition, some regions will
share similar parameters regardless of their geographical distances, due to
the similarities of regions’ demographical information such as income
distribution (Ma et al., 2020; Hu et al., 2020), food environment index, air
pollution (Zhao et al., 2020), and etc.. A spatially varying pattern for
${{\mbox{\boldmath$\theta$}}}(s_{i}),\,i=1,\ldots,n$ is not always valid.
Based on the homogeneity pattern, we focus on the clustering of spatially-
varying parameters. In our setting, we assume that the $n$ parameter vectors
can be clustered into $k$ groups, i.e.,
${\mbox{\boldmath$\theta$}}(s_{i})={\mbox{\boldmath$\theta$}}_{z_{i}}$ where
$z_{i}\in\\{1,2,\ldots,k\\}$.

###  2.2 Geographically Weighted Chinese Restaurant Process

A latent clustering structure can be introduced to accommodate the spatial
heterogeneity on parameters of sub-areas. Under the frequentist framework, the
clustering problem could be solved in a two-stage approach: first obtain the
estimate of number of clusters, $\widehat{k}$, then detect the optimal
clustering assignment among all possible clusterings of $n$ elements into
$\widehat{k}$ clusters. However, in this approach, the performance of the
estimation of cluster assignments highly relies on the estimated number of
clusters, it may ignore uncertainty in the first stage and cause redundant
cluster assignments. Bayesian nonparametric method is a natural remedy to
simultaneously estimate the number of clusters and cluster assignments. The
Chinese restaurant process (CRP; Pitman, 1995; Neal, 2000) offers choices to
allow uncertainty in the number of clusters by assigning a prior distribution
on $(z_{1},z_{2},\ldots,z_{n})$. In CRP, $z_{i},~{}i=2,\ldots,n$ are defined
through the following conditional distribution (also called a Pólya urn
scheme, Blackwell et al., 1973).

$\displaystyle P(z_{i}=c\mid
z_{1},\ldots,z_{i-1})\propto\begin{cases}\absolutevalue{c},&\text{at an
existing cluster labeled}\,c,\\\ \alpha,&\text{at a new cluster}.\end{cases}$
(6)

Here $\absolutevalue{c}$ refers to the size of cluster labeled $c$, and
$\alpha$ is the concentration parameter of the underlying Dirichlet process.
Based on the Pólya urn scheme shown in (6), the customers will have no
preference for sitting with different customers. For the spatial survival
data, nearby regions will share similar environmental effects such as P.M.
2.5, water quality, etc.. These similar effects will lead the nearby sub-
regions to share similar parameters. In order to consider similar effects
caused by geographical distance, we modify the traditional CRP to
geographically weighted CRP (gwCRP) so that the customer will have higher
probability sitting with their familiar customers which are geographically
nearby. We have the conditional distribution of $\bm{\theta}(s_{i})$ given
${\mbox{\boldmath$\theta$}}(s_{1}),\ldots,{\mbox{\boldmath$\theta$}}(s_{i-1})$
based on following definition.

######  Definition 1.

If $G_{0}$ is a continuous distribution and $i>1$, the distribution of
$\bm{\theta}(s_{i})$ given $\bm{\theta}(s_{1}),\ldots,\bm{\theta}(s_{i-1})$ is
proportional to

$\displaystyle
f(\bm{\theta}(s_{i})\mid{\mbox{\boldmath$\theta$}}(s_{1}),\ldots,{\mbox{\boldmath$\theta$}}(s_{i-1}))\propto\sum_{r=1}^{K^{*}}\sum_{j=1}^{i-1}w_{ij}\bm{1}(\bm{\theta}(s_{j})=\bm{\theta}^{*}_{r})\delta_{\bm{\theta}^{*}_{r}}(\bm{\theta}(s_{i}))+\alpha
G_{0}(\bm{\theta}(s_{i})),$ (7)

where $f(\cdot)$ is the distribution density function, $K^{*}$ denote the
number of clusters excluding the $i$-th observation,
$\bm{\theta}^{*}_{1},\ldots,\bm{\theta}^{*}_{K^{*}}$ are $K^{*}$ distinguished
values of $\bm{\theta}_{1},\ldots,\bm{\theta}_{i-1}$, $w_{ij}$ is geographical
weight which is calculated by the distance between $s_{i}$ and $s_{j}$, and
$\delta(\cdot)$ is the Dirac measure.

Based on the Definition 1, we have similar Pólya urn scheme called gwCRP for
conditional distribution in (7) with CRP.

######  Proposition 2.

A Pólya urn scheme of gwCRP is defined as

$\displaystyle P(z_{i}=c\mid
z_{1},\ldots,z_{i-1})\propto\begin{cases}\absolutevalue{c^{*}},&\text{at an
existing cluster labeled}\,c,\\\ \alpha,&\text{at a new cluster}.\end{cases}$
(8)

$\absolutevalue{c^{*}}=\sum_{j=1}^{i-1}w_{ij}{\bf{{1}}}(z_{j}=c)$, where
$w_{ij}$ is the geographical weight.

Compared the existing geographically weighted regression literatures, our
weights are obtained by graph distance between different areas. Following Xue
et al. (2020), we denote a graph as $G$, with set of vertices
$V(G)=\\{v_{1},\ldots,v_{n}\\}$, and set of edges
$E(G)=\\{e_{1},\ldots,e_{m}\\}$. The graph distance between two vertices
$v_{i}$ and $v_{j}$ is defined as:

$d_{v_{i}v_{j}}=\begin{cases}|V(e)|,&\text{if }e\text{ is the shortest path
connecting }v_{i}\text{ and }v_{j},\\\ \infty,&\text{if }v_{i}\text{ and
}\,v_{j}\text{ are not connected},\end{cases}$ (9)

where $|V(e)|$ represents the cardinality of edges in $e$. For the county
level data, we construct the graph $G$ based on adjacency matrix among
different counties. We treat $n$ counties as $n$ vertices of this graph and
$v_{i}$ and $v_{j}$ are connected when the corresponding counties share the
boundary. Based on the graph distance calculated by (9), we calculate the
geographical weights by:

$w_{ij}=\begin{cases}1,&\text{ if }d_{v_{i}v_{j}}\leq 1,\\\
\exp(-d_{v_{i}v_{j}}\times h),&\text{ if }1<d_{v_{i}v_{j}},\end{cases}$ (10)

where $d_{v_{i}v_{j}}$ is the graph distance between areas $i$ and $j$. For
the weighting function in (10), we give the largest weight ($w_{ij}\equiv 1$)
for the areas sharing the same boundaries, which follows the first law of
geography (Tobler, 1970). For simplicity, we refer to gwCRP introduced above
as $\text{gwCRP}(\alpha,h)$, where $\alpha$ is the concentration parameter for
Dirichlet distribution and $h$ is the spatial smoothness parameter.

######  Remark 1.

Based on the Pólya urn scheme defined in (8) and geographical weighting scheme
defined in (10), we find that (i) when $h=0$, the gwCRP reduces to traditional
CRP, which leads to over-clustering problem in estimating of the number of
clusters; (ii) when $h\rightarrow\infty$, a new customer just only choose the
table representing spatially contiguous regions. This will also lead to the
same over-clustering problem as CRP.

###  2.3 gwCRP for Piecewise Constant Hazards Models

Adapting gwCRP to the piecewise constant hazards models, our model and prior
can be expressed hierarchically as:

$\displaystyle\begin{split}&\log\mathcal{L}({\mbox{\boldmath$\beta$}}_{z_{i}},{\mbox{\boldmath$\lambda$}}_{z_{i}},i=1,\ldots,n\mid\bm{D})\\\
&=\sum_{i=1}^{n}\left\\{\sum_{j=1}^{J}d_{ji}\log\lambda_{jz_{i}}+\sum_{\ell=1}^{n_{i}}\delta_{\ell
i}X_{\ell}(s_{i})^{\top}{\mbox{\boldmath$\beta$}}_{z_{i}}-\sum_{j=1}^{J}\lambda_{jz_{i}}\left[\sum_{\ell=1}^{n_{i}}\Delta_{j}(T_{\ell
i})\exp(X_{\ell}(s_{i})^{\top}{\mbox{\boldmath$\beta$}}_{z_{i}})\right]\right\\},\\\
&z_{i}\mid\bm{\pi},k\sim\text{Multinomial}(\pi_{1},\cdots,\pi_{k}),\\\
&\bm{\pi}\sim\text{gwCRP}(\alpha,h),\\\
&\bm{\theta}_{r}\sim\mbox{MVN}(0,\Sigma_{0}),\quad r=1,\ldots,k,\end{split}$
(11)

where
${\mbox{\boldmath$\theta$}}_{r}=(\beta_{1r},\ldots,\beta_{pr},\log\lambda_{1r},\ldots,\log\lambda_{Jr})^{\top}$
is a $p+J$ dimensional vector. And let $k\rightarrow\infty$, and $\Sigma_{0}$
be hyperparameter for base distribution of
${\mbox{\boldmath$\theta$}}_{1},\ldots,{\mbox{\boldmath$\theta$}}_{r}$. We
choose $\Sigma_{0}=100\bm{I}$ in all the simulations and real data analysis
providing noninformative priors. The concentration parameter $\alpha$ controls
the probability of introducing a new cluster which is similar with CRP.
Different values of $h$ lead to different weighting scale for different sub-
regions. In our following simulations and real data analysis, we fix
$\alpha=1$ and tune $h$ with different values.

##  3 Bayesian Inference

In this section, we will introduce the MCMC sampling algorithm, post MCMC
inference method, and Bayesian model selection criterion.

Our goal is to sample from the posterior distribution of the unknown
parameters $k$, $\bm{z}=(z_{1},...,z_{n})\in\\{1,...,k\\}^{n}$,
$\bm{\beta}=(\bm{\beta}_{1},\ldots,\bm{\beta}_{k})$, and
$\bm{\lambda}=(\bm{\lambda}_{1},\ldots,\bm{\lambda}_{k})$. Based on
Proposition 1 and Proposition 2, we can efficiently cycles through the full
conditional distributions of $z_{i}|z_{-i}$ for $i=1,2,\ldots,n$ and
$\bm{\beta}^{\top},\log\bm{\lambda}^{\top}$, where
$z_{-i}=\bm{z}\setminus{\\{z_{i}}\\}$. The marginalization over $k$ can avoid
complicated reversible jump MCMC algorithms or even allocation samplers. The
full conditionals of $z_{1},\ldots,z_{n}$ are given in Proposition 3. The
details of sampling algorithm is given in Section C of Supporting information.

######  Proposition 3.

The full conditional distributions $P(z_{i}=c\mid
z_{-i},\widehat{\bm{\theta}},{\mbox{\boldmath$\theta$}})$ of
$z_{1},\ldots,z_{n}$ is given as

$\propto\left\\{\begin{array}[]{ll}&\left(\sum_{j\neq
i}w_{ij}{\bf{{1}}}(z_{j}=c)\right)(2\pi)^{-\frac{p}{2}}|\widehat{\Sigma}_{i}|^{-\frac{1}{2}}\exp\left\\{-\frac{1}{2}\left((\widehat{\bm{\theta}}(s_{i})-{\mbox{\boldmath$\theta$}}_{c})^{\top}\widehat{\Sigma}_{i}^{-1}(\widehat{\bm{\theta}}(s_{i})-{\mbox{\boldmath$\theta$}}_{c})\right)\right\\}\text{
at existing $c$}\\\
&\alpha(2\pi)^{-\frac{p}{2}}|\widehat{\Sigma}_{i}|^{-\frac{1}{2}}|\Sigma_{0}|^{-\frac{1}{2}}|\widehat{\Sigma}_{i}^{-1}+\Sigma_{0}^{-1}|^{-\frac{1}{2}}\exp\left\\{-\frac{1}{2}\left(\widehat{\bm{\theta}}(s_{i})^{\top}(\widehat{\Sigma}_{i}+\Sigma_{0})^{-1}\widehat{\bm{\theta}}(s_{i})\right)\right\\}\text{
if $c$ is a new cluster}\\\ \end{array}\right.$

where $\widehat{\Sigma}_{i}$ is the estimated variance-covariance matrix of
MLE $\widehat{\bm{\theta}}(s_{i})$ for $i=1,\ldots,n$, and $\Sigma_{0}$ is the
variance hyperparameter for the base distribution of
${\mbox{\boldmath$\theta$}}_{1},\ldots,{\mbox{\boldmath$\theta$}}_{r}$.

We carry out posterior inference on the group memberships $z_{1},\ldots,z_{n}$
by using Dahl’s method (Dahl, 2006), which proceeds as follows

  1. 1.

Define membership matrices
$B^{(l)}=(B(i,j))_{i,j\in\left\\{1,\ldots,n\right\\}}=(\bm{1}(z_{i}^{(l)}=z_{j}^{(l)}))_{n\times
n}$, where $l=1,\ldots,B$ indexes the number of retained MCMC draws after
burn-in, and $\bm{1}(\cdot)$ is the indicator function.

  2. 2.

Calculate the average membership matrix
$\overline{B}=\frac{1}{B}\sum_{l=1}^{L}B^{(l)}$, where the summation is
element-wise.

  3. 3.

Identify the most _representative_ posterior sample as the one that is closest
to $\overline{B}$ with respect to the element-wise Euclidean distance
$\sum_{i=1}^{n}\sum_{j=1}^{n}(B^{(l)}(i,j)-\overline{B}(i,j))^{2}$ among the
retained $l=1,\ldots,B$ posterior samples.

Therefore, the posterior estimates of cluster memberships $z_{1},\ldots,z_{n}$
and model parameters $\bm{\theta}$ can be obtained based on the draw
identified by Dahl’s method.

We recast the choice of decaying parameter $h$ as a model selection problem.
We use the Logarithm of the Pseudo-Marginal Likelihood (LPML; Ibrahim et al.,
2001) based on conditional predictive ordinate (CPO; Gelfand et al., 1992;
Geisser, 1993; Gelfand and Dey, 1994) to select $h$. The LPML is defined as

$\text{LPML}=\sum_{i=1}^{N}\text{log}(\text{CPO}_{i}),$ (12)

where $\text{CPO}_{i}$ is the i-th conditional predictive ordinate. Following
Chen et al. (2000), a Monte Carlo estimate of the CPO, within the Bayesian
framework, can be obtained as

$\widehat{\text{CPO}}_{i}^{-1}=\frac{1}{B}\sum_{b=1}^{B}\frac{1}{f(D_{i}|{\mbox{\boldmath$\theta$}}_{z_{i}}^{b})},$
(13)

where $B$ is the total number of Monte Carlo iterations,
${\mbox{\boldmath$\theta$}}_{z_{i}}^{b}$ is the $b$-th posterior sample, and
$f(\cdot)$ is the likelihood function define in (3). An estimate of the LPML
can subsequently be calculated as:

$\widehat{\text{LPML}}=\sum_{i=1}^{N}\text{log}(\widehat{\text{CPO}}_{i}).$
(14)

A model with a larger LPML value is preferred.

##  4 Simulation

###  4.1 Simulation Setting and Evaluation Metrics

In this section, we present simulation studies under four different designs to
illustrate the performance of our proposed gwCRP method and compare with
traditional CRP, in terms of both clustering configuration and estimation of
regression coefficients and piecewise constant baseline hazards under
proportional hazards model. Survival datasets that resemble the SEER
respiratory cancer data for Louisiana are generated. The censoring rate is
around 30%. We design four different geographical clustering patterns in
Louisiana state, which are shown in Figure 1. Designs I and III have three
true clusters, and Designs II and IV have two true clusters. In addition,
Designs II and III both have one cluster consisting of two disjoint areas
since, in practice, it is still possible for two distant counties to belong to
the same cluster. Design IV has two clusters both consisting of disjoint
areas.

Figure 1: Geographical clustering patterns in Louisiana state of simulation
designs (This figure appears in color in the electronic version of this
article, and any mention of color refers to that version.)

For each design, 100 replicate datasets are generated under proportional
hazards model with piecewise constant baseline hazard. In each replicate, we
generate survival data of 60 subjects for each county, including three
regression covariates from $N(0,1)$ i.i.d., survival time and censoring. We
set three pieces for the baseline hazards with cutting points 1.5 and 6 for
all designs, and over four designs, we have three true clusters at maximum,
and the true regression coefficients and baseline hazards used are chosen from
${\mbox{\boldmath$\beta$}}_{1}=(1,0.5,1),{\mbox{\boldmath$\lambda$}}_{1}=(0.045,0.036,0.045)$,
${\mbox{\boldmath$\beta$}}_{2}=(1.5,1,1),{\mbox{\boldmath$\lambda$}}_{2}=(0.045,0.036,0.036)$,
and
${\mbox{\boldmath$\beta$}}_{3}=(2,0.5,1.5),{\mbox{\boldmath$\lambda$}}_{3}=(0.036,0.045,0.0495)$.
Censoring times are generated independently by taking the minimum of 150 and
random values from Exp(0.01) with expectation 100. For each replicate, we set
$\alpha=1$ and run different values of $h$, from 0 to 2 with grid 0.2, and
from 3 to 10 with grid 1, and select the optimal $h$ via LPML. A total of 2000
MCMC iterations are run for each replicate, with the first 500 iterations as
burn-in.

To compare the performance of clustering of gwCRP under different values of
$h$, both estimation of the number of clusters and the matchability of
clustering configurations are reported. In our simulation, we use mean Rand
Index (Rand, 1971) which is obtained by using R-package fossil (Vavrek, 2011)
to measure the clustering performance.

In addition to clustering performance, we further evaluate the estimation
performance of covariates coefficients and baseline hazards, which is assessed
by average of bias (AB) and average of mean squared error (AMSE) defined as
follows. Let $\bm{z}=(z_{1},\ldots,z_{n})$ be the true clustering label
vector, ${\mbox{\boldmath$\theta$}}_{r}(s_{i})$ be the true parameter value of
cluster $r$, $\kappa_{r}=\sum_{i=1}^{64}{\bf{{1}}}(z_{i}=r)$ be the number of
counties in cluster $r$, $r=1,\ldots,k$, $\sum_{r=1}^{k}\kappa_{r}=n$, and for
simulated data set $t$, let
$\widehat{{\mbox{\boldmath$\theta$}}}_{(t)}(s_{i})$ be Dahl’s method estimate
at location $s_{i}$ for $t$-th replicates. Then AB is calculated as

$\text{AB}=\frac{1}{k}\sum_{r=1}^{k}\frac{1}{\kappa_{r}}\sum_{i|z_{i}=r}\frac{1}{100}\sum_{t=1}^{100}(\widehat{{\mbox{\boldmath$\theta$}}}_{(t)}(s_{i})-{\mbox{\boldmath$\theta$}}_{r}(s_{i})),$

and AMSE is calculated as

$\text{AMSE}=\frac{1}{k}\sum_{r=1}^{k}\frac{1}{\kappa_{r}}\sum_{i|z_{i}=r}\frac{1}{100}\sum_{t=1}^{100}(\widehat{{\mbox{\boldmath$\theta$}}}_{(t)}(s_{i})-{\mbox{\boldmath$\theta$}}_{r}(s_{i}))^{2},$

which calculates mean squared errors for each cluster first, and then average
across clusters.

###  4.2 Simulation Results

Figure 2 shows the histogram of $k$ estimates and boxplots of Rand Index under
different $h$ and the optimal selected by LPML for four simulation designs. We
see that when $h=0$, the proposed gwCRP method is identical to the traditional
CRP method, and in this case, CRP always tends to over-cluster and often
yields smaller Rand Index than the results under $h>0$. Another important
trend is that, as $h$ increases, the estimated number of clusters decreases
first and then increases, and the Rand Index increases first and then
decreases as $h$ becomes too large. As we discussed in Remark 1, this is
because when $h$ increases from 0, the spatial patterns in the data is
captured by the proposed gwCRP method. However, as $h\rightarrow\infty$, the
geographical weights $w_{ij}$ for spatial-discontiguous counties become 0,
which means only adjacent counties can be classified into the same cluster,
therefore leading to over-clustering phenomenon again. It is also discovered
that the clustering perfomance under optimal $h$ selected by LPML is very
well, with the probability of selecting true number of clusters always greater
than 0.75, and Rand Index larger than or similar to the highest results
attained by some fixed value of $h$.

Figure 2: Histogram of estimates of $k$ and boxplot of Rand Index under
different $h$ and LPML selection for simulation designs. The average h
selected by LPML is 1.296 in Design 1, 1.412 in Design 2, 1.366 in Design 3,
1.602 in Design 4.

Table 1 summarizes the AB and AMSE results of estimating parameters of gwCRP
under different $h$ for different designs. For simplicity of summary results,
here the AMSE of $\beta$ is the average of AMSE of
$\beta_{1},\beta_{2},\beta_{3}$ since they have similar scales, and the value
of $\log{\mbox{\boldmath$\lambda$}}$ is the average of AMSE of
$\log\lambda_{1},\log\lambda_{2},\log\lambda_{3}$, respectively.

For Designs II, III and IV, the absolute value of AB decrease as $h$ increase
from 0 to moderate values, and increase again as $h$ increase to relatively
large values. For all four designs, the absolute values of AB for $\lambda$’s
of optimal $h$ selected by LPML always are the smallest, and the absolute
values for both $\beta$’s and $\lambda$’s of optimal $h$ selected by LPML are
always smaller than the values of traditional CRP. The patterns in AMSE are
more clear when comparing different methods, that traditional CRP has the
largest AMSE and AMSE decrease as $h$ increase from 0 to moderate values, and
increase again as $h$ increase to relatively large values. The results of
optimal $h$ selected by LPML also has the best performance in estimation.

Table 1: AB and AMSE for parameter estimation under different $h$ and LPML selection for different simulation designs Method | Parameter | Design I | Design II | Design III | Design IV
---|---|---|---|---|---
|  | AB | AMSE | AB | AMSE | AB | AMSE | AB | AMSE
gwCRP $h=0.6$ | $\beta$ | -0.0063 | 0.0069 | 0.0038 | 0.0067 | -0.0065 | 0.0083 | -0.0027 | 0.0060
| $\lambda$ | 0.0815 | 0.0193 | 0.0732 | 0.0186 | 0.0797 | 0.0219 | 0.0827 | 0.0190
gwCRP $h=1.2$ | $\beta$ | -0.0078 | 0.0065 | -0.0002 | 0.0058 | -0.0061 | 0.0087 | -0.0017 | 0.0049
| $\lambda$ | 0.0818 | 0.0194 | 0.0775 | 0.0171 | 0.0794 | 0.0217 | 0.0815 | 0.0172
gwCRP $h=2$ | $\beta$ | -0.0096 | 0.0068 | -0.0006 | 0.0055 | -0.0056 | 0.0085 | -0.0051 | 0.0049
| $\lambda$ | 0.0849 | 0.0190 | 0.0814 | 0.0158 | 0.0798 | 0.0216 | 0.0873 | 0.0191
gwCRP $h=6$ | $\beta$ | -0.0061 | 0.0129 | 0.0030 | 0.0072 | -0.0299 | 0.0204 | -0.0064 | 0.0074
| $\lambda$ | 0.0805 | 0.0281 | 0.0770 | 0.0217 | 0.1005 | 0.0296 | 0.0852 | 0.0195
gwCRP Optimal | $\beta$ | -0.0039 | 0.0059 | 0.0042 | 0.0055 | 0.0074 | 0.0067 | -0.0005 | 0.0035
| $\lambda$ | 0.0661 | 0.0177 | 0.0711 | 0.0145 | 0.0732 | 0.0203 | 0.0777 | 0.0177
CRP | $\beta$ | -0.0046 | 0.0086 | 0.0018 | 0.0092 | -0.0056 | 0.0089 | 0.0003 | 0.0082
| $\lambda$ | 0.0760 | 0.0228 | 0.0717 | 0.0233 | 0.0787 | 0.0239 | 0.0742 | 0.0223

A sensitivity analysis regarding $\alpha$ and the weighting function is
conducted. $\alpha=0.5,1,2,5$ and the weighting function
$w_{ij}=\exp(-d_{v_{i}v_{j}}^{2}\times
h^{2}){\bf{{1}}}\\{d_{v_{i}v_{j}}>1\\}+{\bf{{1}}}\\{d_{v_{i}v_{j}}\leq 1\\}$
which has a faster decay to 0 are ran, and all the results are presented in
Section D of Supporting information. The results show that the results of
optimal $h$ selected by LPML are insensitive to the choice of $\alpha$ and
weighting function.

In a brief conclusion based on our simulation studies, gwCRP models have
better performance than CRP both for clustering and parameter estimation. Our
proposed model selection criterion, LPML, can nearly select the best
performing $h$ value for both clustering and parameter estimation.

##  5 SEER Respiratory Cancer Data

###  5.1 Data Description

In this section, we apply our proposed model to analyze respiratory cancer
data in Louisiana state, which is downloaded from the Surveillance,
Epidemiology, and End Results (SEER) Program.We analyzed the survival time of
respiratory cancer patients using the SEER public use data (SEER 1973-2016
Public-Use). We refer to Mu et al. (2020) for the detailed data clean
description. After cleaning, there are 16213 observations left, and the
censoring rate is 30.44%. We select Age, Gender, Cancer grade and Historical
stage of cancer for our analysis, and give the summary of survival times and
covariates in Table 2. The median survival times for patients in each county
are provided in Section E of Supporting information .

Table 2: Demographics for the studied dataset. For continuous variables, the mean and standard deviation (SD) are reported. For binary variables, the frequency and percentage of each class are reported. | Mean(SD) / Frequency (Percentage)
---|---
Survival Time | 22.43 (31.90)
Event | 12.63 (18.32)
Censor | 44.85 (43.06)
Diagnostic Age | 66.55 (11.66)
Sex |
Female | 6548 ($40.39\%$)
Male | 9665 ($59.61\%$)
Cancer Grade |
the class of lower grades | 5307 ($32.73\%$)
the class of III or IV | 10906 ($67.27\%$)
Historical Stage 111Distant stage means that a tumor has spread to areas of the body distant or remote from the primary tumor |
not distant | 9005 ($55.54\%$)
distant | 7208 ($44.46\%$)

We first fit the Cox model of patients for each county using the covariates
selected. The regression coefficients are visualized in Section E of
Supporting information. From results shown in Supporting information, it is
seen that some counties have similar characteristics, no limited to only
adjacent counties, indicating possibilities of globally discontiguous
clusters.

###  5.2 Data Analysis

To select the optimal number of pieces for the baseline hazard and $h$, we run
$J=2$ with the cutpoints $(0,9.01)$, $J=3$ with cutpoints $(0,3.01,9.01)$,
$J=4$ with cutpoints $(0,1.01,4.01,9.01)$, and $J=5$ with cutpoints
$(0,1.01,3.01,5.01,9.01)$. The cutpoints are set by dividing the start point
$0$ and the median survival time $9.01$ by quantiles evenly to ensure there
are events at each piece for each county. For each $J$, we run $h$ from 0 to
10 with grid 0.1, and for each combination of $J$ and $h$, 5000 MCMC
iterations are run and drop the first 2000 as burn-in. The optimal values
selected by LPML is $J=4$ and $h=9.0$, under which the corresponding estimate
of number of clusters is two, while the traditional CRP classifies the
counties into five clusters. The trace plots of different chains of posterior
samples of the estimates for selected counties are presented in the Section E
of Supporting information to show the convergence of the MCMC. The plots of
clustering patterns of CRP and gwCRP Optimal are shown in Figure 3, from which
it is seen that the gwCRP captures the globally discontiguous clusters very
well. The estimates and $95\%$ Credible Intervals of regression covariates
coefficients and baseline hazards obtained by gwCRP Optimal are given in Table
3, from which we see that, though the baseline hazards are similar, the
regression covariates coefficients are quite different across different
clusters. We see that our proposed method successfully detects both spatially
contiguous cluster and discontinuous cluster simultaneously. The parameter
estimates for Age are positive in all counties, indicating that older patients
on average are more likely to have the event than younger patients. For the
counties in cluster 1, the diagnostic ages has higher hazards effects than
other counties. However, for the counties in Cluster 2, male, later cancer
stage will have higher hazards effects than other counties. The historical
distance stage effects are very similar in two clusters which indicates that
the subjects with tumor spreading will have higher hazards.

Figure 3: Clustering patterns of counties in Louisiana state under CRP (when J=4) and gwCRP Optimal (when J=4, h=9.0) methods (This figure appears in color in the electronic version of this article, and any mention of color refers to that version.) Table 3:  Estimates results of regression coefficients and baseline hazards obtained by gwCRP Optimal ($J=4,~{}h=9.0$). The $95\%$ Credible Interval for estimates of Cluster One is calculated by the $95\%$ HPD Interval of County 13, and The $95\%$ Credible Interval for estimates of Cluster Two is calculated by the $95\%$ HPD Interval of County 33, where the counties were selected by the minimum Euclidean distance from the posterior mean to the average estimate.  Parameter | Cluster 1 | Cluster 2
---|---|---
Estimate |  $95\%$ Credible Interval | Estimate |  $95\%$ Credible Interval
$\beta_{\text{Age}}$ | 0.1847 | (0.1693, 0.2158) | 0.0728 | (0.0138, 0.2056)
$\beta_{\text{Sex}}$ | 0.1239 | (0.0912, 0.1732) | 0.3411 | (0.1118, 0.5189)
$\beta_{\text{Grade}}$ | 0.5075 | (0.4366, 0.5291) | 0.8290 | (0.4583, 0.9578)
$\beta_{\text{Hist-Stage}}$ | 1.3271 | (1.2926, 1.3824) | 1.4434 | (1.3150, 1.7101)
$\lambda_{1}$ | 1.0690 | (0.9938, 1.0912) | 1.0359 | (0.9060, 1.2330)
$\lambda_{2}$ | 1.0716 | (0.9909, 1.1172) | 1.0877 | (0.8432, 1.3035)
$\lambda_{3}$ | 0.9843 | (0.9787, 1.1040) | 1.0912 | (0.8535, 1.2899)
$\lambda_{4}$ | 1.0040 | (0.9583, 1.0513) | 1.0209 | (0.8724, 1.2091)

##  6 Discussion

In this paper, we proposed a geographically weighted Chinese restaurant
process to capture spatial homogeneity of regression coefficients and baseline
hazards based on piecewise constant hazard model. An efficient MCMC algorithm
is proposed for our methods without complicated reversible jump algorithm.
Extensive simulation results are carried out to show that our proposed method
has better clustering performance than the traditional CRP in spatial
homogeneity pursuit for survival data. Simulation studies also show that our
proposed methods have promising results in coefficients and baseline hazard
estimation. An application to analysis of SEER data provides an interesting
illustration of our proposed methods.

Furthermore, four topics beyond the scope of this paper are worth further
investigation. In this paper, our proposed algorithm is based on two step
estimation under piecewise constant proportional hazard model assumption.
Proposing an efficient sampling algorithm without Laplace approximation is an
important future work. Furthermore, we fixed the number of pieces of baseline
hazards in both simulation studies and real data analysis. Imposing adaptive
number of pieces model in baseline hazards is devoted for future research. In
addition, variable selection approaches based on hierarchical CRP (Griffiths
et al., 2004) is also worth being investigated. Allowing different covariates
and baseline hazard share different clustering processes is also an important
future work.

## Acknowledgement

The authors would like to thank the editor, the associate editors and two
reviewers for their valuable comments which help improve the presentation of
this paper.

## Data Availability Statement

The data that support the findings of this paper are available from the
corresponding author upon reasonable request.

## References

  * Banerjee et al. (2014) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2014).  Hierarchical modeling and analysis for spatial data.  Crc Press.
  * Banerjee and Dey (2005) Banerjee, S. and Dey, D. K. (2005).  Semiparametric proportional odds models for spatially correlated survival data.  Lifetime Data Analysis 11, 175–191.
  * Banerjee et al. (2003) Banerjee, S., Wall, M. M., and Carlin, B. P. (2003).  Frailty modeling for spatially correlated survival data, with application to infant mortality in Minnesota.  Biostatistics 4, 123–142.
  * Bhatt and Tiwari (2014) Bhatt, V. and Tiwari, N. (2014).  A spatial scan statistic for survival data based on weibull distribution.  Statistics in medicine 33, 1867–1876.
  * Blackwell et al. (1973) Blackwell, D., MacQueen, J. B., et al. (1973).  Ferguson distributions via pólya urn schemes.  The Annals of Statistics 1, 353–355.
  * Blei and Frazier (2011) Blei, D. M. and Frazier, P. I. (2011).  Distance dependent Chinese restaurant processes.  Journal of Machine Learning Research 12, 2461–2488.
  * Chen et al. (2000) Chen, M.-H., Qi-Man, S., and G, I. J. (2000).  Monte Carlo Methods in Bayesian Computation.  New York: Springer-Verlag.
  * Cox (1972) Cox, D. R. (1972).  Regression models and life-tables.  Journal of the Royal Statistical Society B 34, 187–220.
  * Dahl (2006) Dahl, D. B. (2006).  Model-based clustering for expression data via a Dirichlet process mixture model.  Bayesian inference for gene expression and proteomics 4, 201–218.
  * Friedman et al. (1982) Friedman, M. et al. (1982).  Piecewise exponential models for survival data with covariates.  The Annals of Statistics 10, 101–113.
  * Geisser (1993) Geisser, S. (1993).  Predictive Inference: An Introduction.  London: Chapman & Hall.
  * Gelfand and Dey (1994) Gelfand, A. E. and Dey, D. K. (1994).  Bayesian model choice: asymptotics and exact calculations.  Journal of the Royal Statistical Society: Series B (Methodological) 56, 501–514.
  * Gelfand et al. (1992) Gelfand, A. E., Dey, D. K., and Chang, H. (1992).  Model determination using predictive distributions with implementation via sampling-based-methods (with discussion).  In In Bayesian Statistics 4. University Press.
  * Gelfand et al. (2003) Gelfand, A. E., Kim, H.-J., Sirmans, C., and Banerjee, S. (2003).  Spatial modeling with spatially varying coefficient processes.  Journal of the American Statistical Association 98, 387–396.
  * Green (1995) Green, P. J. (1995).  Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  Biometrika 82, 711–732.
  * Griffiths et al. (2004) Griffiths, T. L., Jordan, M. I., Tenenbaum, J. B., and Blei, D. M. (2004).  Hierarchical topic models and the nested Chinese restaurant process.  In Advances in neural information processing systems, pages 17–24.
  * Henderson et al. (2012) Henderson, R., Shimakura, S., and Gorst, D. (2012).  Modeling spatial variation in leukemia survival data.  Journal of the American Statistical Association .
  * Hu et al. (2020) Hu, G., Geng, J., Xue, Y., and Sang, H. (2020).  Bayesian spatial homogeneity pursuit of functional data: an application to the U.S. income distribution.
  * Hu and Huffer (2020) Hu, G. and Huffer, F. (2020).  Modified Kaplan–Meier estimator and Nelson–Aalen estimator with geographical weighting for survival data.  Geographical Analysis 52, 28–48.
  * Hu et al. (2020) Hu, G., Xue, Y., and Huffer, F. (2020).  A comparison of Bayesian accelerated failure time models with spatially varying coefficients.  Sankhya B pages 1–17.
  * Huang et al. (2007) Huang, L., Pickle, L. W., Stinchcomb, D., and Feuer, E. J. (2007).  Detection of spatial clusters: application to cancer survival as a continuous outcome.  Epidemiology 18, 73–87.
  * Ibrahim et al. (2001) Ibrahim, J. G., Chen, M.-H., and Sinha, D. (2001).  Bayesian Survival Analysis.  New York: Springer-Verlag.
  * Lee et al. (2017) Lee, J., Gangnon, R. E., and Zhu, J. (2017).  Cluster detection of spatial regression coefficients.  Statistics in Medicine 36, 1118–1133.
  * Lee et al. (2019) Lee, J., Sun, Y., and Chang, H. H. (2019).  Spatial cluster detection of regression coefficients in a mixed-effects model.  Environmetrics page e2578.
  * Li and Sang (2019) Li, F. and Sang, H. (2019).  Spatial homogeneity pursuit of regression coefficients for large datasets.  Journal of the American Statistical Association pages 1–21.
  * Lu et al. (2018) Lu, J., Li, M., and Dunson, D. (2018).  Reducing over-clustering via the powered Chinese restaurant process.  arXiv preprint arXiv:1802.05392 .
  * Ma et al. (2020) Ma, Z., Xue, Y., and Hu, G. (2020).  Heterogeneous regression models for clusters of spatial dependent data.  Spatial Economic Analysis 15, 459–475.
  * Miller and Harrison (2013) Miller, J. W. and Harrison, M. T. (2013).  A simple example of Dirichlet process mixture inconsistency for the number of components.  In Advances in Neural Information Processing Systems, pages 199–206.
  * Mu et al. (2020) Mu, J., Liu, Q., Kuo, L., and Hu, G. (2020).  Bayesian variable selection for cox regression model with spatially varying coefficients with applications to louisiana respiratory cancer data.  arXiv preprint arXiv:2008.00615 .
  * Neal (2000) Neal, R. M. (2000).  Markov chain sampling methods for Dirichlet process mixture models.  Journal of Computational and Graphical Statistics 9, 249–265.
  * Pitman (1995) Pitman, J. (1995).  Exchangeable and partially exchangeable random partitions.  Probability Theory and Related Fields 102, 145–158.
  * Rand (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association 66, 846–850.
  * SEER (2016) SEER, P. (2016).  Public-use data (1973-2015). national cancer institute, dccps, surveillance research program, cancer statistics branch, released april 2016, based on the november 2015 submission.
  * Tobler (1970) Tobler, W. R. (1970).  A computer movie simulating rrban growth in the Detroit region.  Economic Geography 46, 234–240.
  * Vavrek (2011) Vavrek, M. J. (2011).  Fossil: palaeoecological and palaeogeographical analysis tools.  Palaeontologia Electronica 14, 16\.
  * Xue et al. (2020) Xue, Y., Schifano, E. D., and Hu, G. (2020).  Geographically weighted Cox regression for prostate cancer survival data in Louisiana.  Geographical Analysis 52, 570–587.
  * Zhang and Lawson (2011) Zhang, J. and Lawson, A. B. (2011).  Bayesian parametric accelerated failure time spatial model and its application to prostate cancer.  Journal of Applied Statistics 38, 591–603.
  * Zhao et al. (2020) Zhao, P., Yang, H.-C., Dey, D. K., and Hu, G. (2020).  Bayesian spatial homogeneity pursuit regression for count value data.  arXiv preprint arXiv:2002.06678 .
  * Zhou et al. (2008) Zhou, H., Lawson, A. B., Hebert, J. R., Slate, E. H., and Hill, E. G. (2008).  Joint spatial survival modeling for the age at diagnosis and the vital outcome of prostate cancer.  Statistics in Medicine 27, 3612–3628.

## Supporting Information

Web Appendices, Tables, and Figures referenced in Sections 2-5 and R scripts
for simulationsand real data examples are available with this paper at the
Biometrics website on WileyOnline Library. The R code for the computations of
this paper is available at https://github.com/lj-geng/GWCRP.