Bayesian Spatial Split-Population Survival Model with Applications to Democratic Regime Failure and Civil War Recurrence

Abstract

The underlying risk factors associated with the duration and termination of biological, sociological, economic, or political processes often exhibit spatial clustering. However, existing nonspatial survival models, including those that account for “immune” and “at-risk” subpopulations, assume that these baseline risks are spatially independent, leading to inaccurate inferences in split-population survival settings. In this paper, we develop a Bayesian spatial split-population survival model that addresses these methodological challenges by accounting for spatial autocorrelation among units in terms of their probability of becoming immune and their survival rates. Monte Carlo experiments demonstrate that, unlike nonspatial models, this spatial model provides accurate parameter estimates in the presence of spatial autocorrelation. Applying our spatial model to data from published studies on authoritarian reversals and civil war recurrence reveals that accounting for spatial autocorrelation in split-population models leads to new empirical insights, reflecting the need to theoretically and statistically account for space and non-failure inflation in applied research.

1. Introduction

Originally used to study human survival rates following the onset of a disease or the administration of medical treatment, parametric and semi-parametric tools for modeling time-to-event data or “survival times” have been used to study innumerable biological, industrial, psychological, social, and political phenomena. However, two common types of heterogeneity in the data generation process (d.g.p.) of many time-to-event applications violate the core assumptions of conventional survival models. The first is the presence of non-failure cases resulting from “immunity” to a failure event or being “cured” from that event due to some treatment. Cases that will never experience the event of interest violate the assumption that all right-censored observations eventually experience the failure event even if the failure is not observed. The frequent need to relax this assumption in applied settings has given rise to a class of split-population (SP) survival models that first estimate the probability of being immune or at risk of experiencing the event and subsequently estimate the time until that event occurs, conditional upon not being immune to the event. In other words, SP survival models do not assume that every observation will eventually experience the event. “Instead, the model splits the population into two groups—one that will experience the event and one that will not” [1] (p. 148). The probability of a case being immune to the event is estimated in the first (split) stage as a binary process modeled with a specified set of covariates, then the survival stage is modeled with a specified baseline function representing the time until those cases at risk of experiencing the event actually do so, again conditional upon covariates.

These tools have been useful for modeling a wide range of phenomena, including oncological studies of the survival of breast cancer patients [2] and melanoma relapse [3], the occurrence of interstate war [4], susceptibility to and mortality from parasitic infection among river salmon [5], and criminal recidivism [6]. SP survival models themselves have been extended to incorporate independent and identically distributed (i.i.d.) frailties [7], to account for random right-censoring [8], to address misclassified failures [9], and to account for “triadic duration” independence [10]; however, existing formulations often ignore the effects of spatial clustering among units, at least in the first stage [11].

The spatial clustering of common unobserved characteristics among units that may affect their baseline risk of experiencing a failure event violates a second core assumption of conventional survival models, namely, that units are conditionally independent. Such spatial autocorrelation differs from spatial dependence in that it cannot be accounted for with i.i.d. frailty terms or spatial lags [12,13]. Spatially weighted frailties have been appropriately incorporated into conventional survival models via Bayesian estimation [12,14,15,16]. However, existing spatial survival models cannot accommodate any heterogenous mixture of immune and at-risk populations, nor have they included “covariates and spatial random effects as regressors in the cure rate portion of the model, instead of just the log-relative risk portion” [11] (p. 274).

In this article, we develop a parametric Bayesian spatial split-population (SP) survival model that can incorporate time-varying covariates. Rather than adopting the frequentist maximum likelihood estimation approach to finding parameter values, Bayesian estimation more flexibly utilizes Bayes’ Rule to estimate parameters based on iterated updates to pre-specified priors that define baseline expectations about the probability distribution of the phenomenon of interest. This method is particularly useful for estimating parameters in a split-population model, which can place high demands on the observed data.

Similar to a conventional split-population survival model, our approach consists of a split-stage equation that estimates the probability of a unit being immune from a failure event and a second-stage equation that estimates the survival probability conditioned upon the subject being at risk of failure. In our spatial frailty model, however, each of the two equations may include spatially autocorrelated frailties with a joint distribution that is interpretable in a spatial context. This allows analysts to eschew the assumption that the frailties themselves are i.i.d. The hierarchical model leverages the flexibility of Bayesian estimation using a Markov Chain Monte Carlo (MCMC) sampling algorithm whereby the frailties of “neighboring” units and any spatial autocorrelations among them are incorporated into each equation via a conditionally autoregressive (CAR) prior. MCMC algorithms are tools for estimating analytically complex probability densities by randomly and repeatedly drawing samples from a distribution (Monte Carlo) such that each sample depends on the prior one but not those before it (Markov Chain). Thus, unlike regular SP survival models, the Bayesian spatial SP survival model can account for spatial autocorrelation in a unit’s propensity for being at risk of experiencing an event as well as the time it takes for that event to occur.

After presenting the Bayesian spatial SP survival model and describing the slice-sampling algorithm used for estimation in a Bayesian framework, we illustrate its properties through a series of Monte Carlo experiments. The results reveal that (i) the Bayesian spatial SP model reliably retrieves its true parameter values regardless of the size of the immune fraction or degree of spatial autocorrelation and that (ii) nonspatial models produce biased parameter estimates if the true d.g.p. includes spatial clustering. We then apply the spatial SP survival model to replication data from two prominent studies in political science. The first is a previous application of a nonspatial split-population survival model used to study whether democratic countries consolidate, survive, or revert to a dictatorship [17]. The second is a conventional survival analysis of whether civil wars are more or less likely to recur after terminating [18]. We discuss theoretical reasons to expect spatial clustering in each of these contexts, then demonstrate empirically that spatial clustering does indeed exist in the replication data. Then, we show that applying our Bayesian spatial SP survival model to these data significantly alters the previously reported results. In light of the evidence of spatial autocorrelation in each application, the new results from our spatial model indicate that faulty inferences can result from ignoring spatial heterogeneity when modeling survival processes.

2. (Spatial) Split-Population Survival Model

2.1. Model Development

Supposing that $i=\{1,2,\dots ,N\}$ are the units in the survival data, we define $f\left(t\right)$ as the probability density function and $F\left(t\right)$ as the cumulative distribution function. Thus, $S\left(t\right)=1-F\left(t\right)$ is the survival distribution and $h\left(t\right)=\frac{f\left(t\right)}{S\left(t\right)}$ is the hazard rate. The general likelihood of the conventional survival model is proportional to

$$L=\prod _{i=1}^N{\left[f\left(t_i\right)\right]}^{{\tilde{C}}_i}{\left[S\left(t_i\right)\right]}^{1-{\tilde{C}}_i},$$

(1)

where ${\tilde{C}}_i=1$ are the units that fail and ${\tilde{C}}_i=0)$ are the units that do not fail, and as such are “right-censored.” Two subpopulations can potentially exist in the survival data that researchers use for empirical analysis: an “at-risk” fraction of cases that can fail, and an “immune” fraction of cases that will not fail, implying that units in this fraction do not experience the event of interest [19,20,21]. These two subpopulations are accounted for in split-population survival models (with or without unit-specific frailties) by estimating the probability of a given unit being in the immune fraction and the influence of covariates on the at-risk fraction’s hazard rate [19,22,23].

The split-population survival model for a duration t that splits the sample in the manner described above is constructed as follows. First, we define ${\alpha }_i=Pr(Y_i=1)$ as the probability of units entering the immune fraction, which can be estimated via logit:

$${\alpha }_i=\frac{exp({\mathbf{Z}}_i\gamma +V_i)}{1+exp({\mathbf{Z}}_i\gamma +V_i)}$$

(2)

where ${\mathbf{Z}}_i$ are p2-dimensional covariates, $\gamma$ is the corresponding parameter vector in ${\mathbb{R}}^{p_2}$, and $V_i\sim N(0,{\sigma }^2)$ are the nonspatial unit-specific frailties (random effects), which are assumed to be independent and identically distributed (i.i.d.). The nonspatial i.i.d. unit-specific frailties $V_i$ in the model’s split-stage accounts for unobserved heterogeneity on ${\alpha }_i$ while being independent of other random effects. In the split-population model’s survival stage, however, $W_i\sim N(0,{\sigma }^2)$ denotes the nonspatial i.i.d. unit-specific frailties. As such, these nonspatial frailties capture unobserved factors that potentially influence the units’ distinct risks of experiencing the event of interest. Hence, the proportional hazards of the split-population survival model with nonspatial unit-specific i.i.d. frailties are

$$h\left(t_i\right|{\mathbf{X}}_i\beta ,W_i)=h_0\left(t_i\right){\omega }_{i}exp\left({\mathbf{X}}_i\beta \right)=h_0\left(t_i\right)exp({\mathbf{X}}_i\beta +W_i),$$

(3)

where $h_0\left(t_i\right)$ is the baseline hazard (which can be Weibull, log-logistic, or log-normal distributed), $log{\omega }_i=W_i$, ${\mathbf{X}}_i$ represents the p1-dimensional covariates, and $\beta$ is the corresponding parameter vector in ${\mathbb{R}}^{p_1}$. As discussed below, while we use the Weibull distribution in the presented Monte Carlo experiments and applications, the experiments perform similarly when using either of the other two parametric distributions.

Suppose that we need to incorporate time-varying covariates in our split population. Let $t0$ be the unique “entry time” and let t be the “exit time” for each period. Suppose that j is the beginning of the time period. Then, each unit i’s elapsed time from inception until (i) j can be denoted as $t{0}_{ij}$ and (ii) the end of period j is $t_{ij}$. In this case, ${\tilde{C}}_{ij}=0$ implies that the observation can be censored, while ${\tilde{C}}_{ij}=1$ indicates that the observation has failed at $t_{ij}$. The probability of survival until period j is now $S_i\left(t0\right)=1-F\left(t0\right)$, where $F\left(t0\right)={\int }_0^{t0}f\left(t0\right)$. In this case, both subpopulations contribute to the log-likelihood of the split-population survival model with nonspatial i.i.d. frailties as follows: $S\left(t0\right)=1-F\left(t0\right)$ (where $F\left(t0\right)={\int }_0^{t0}f\left(u\right)du$) is the probability of survival until j. Thus, the log-likelihood of the split-population survival model with nonspatial i.i.d. frailties is proportional to

$$\ln L=\sum _{i=1}^N\left\{{\tilde{C}}_{ij}\ln \left[(1-{\alpha }_{ij})\frac{f\left(t_{ij}\right|{\mathbf{X}}_{ij}\beta ,W_i)}{S\left(t{0}_{ij}\right|{\mathbf{X}}_{ij}\beta ,W_i)}\right]+(1-{\tilde{C}}_{ij})\ln \left[{\alpha }_i+(1-{\alpha }_i)\frac{S\left(t_{ij}\right|{\mathbf{X}}_{ij}\beta ,W_i)}{S\left(t{0}_{ij}\right|{\mathbf{X}}_{ij}\beta ,W_i)}\right]\right\}$$

(4)

where ${\alpha }_{ij}=\frac{exp({\mathbf{Z}}_{ij}\gamma +V_i)}{1+exp({\mathbf{Z}}_{ij}\gamma +V_i)}$ is the split-stage equation; the model’s survival stage estimates the effect of covariates ${\mathbf{X}}_{ij}$ on the probability of survival conditional on each unit being at-risk along with the baseline hazard. Here, $V_i$ and $W_i$ are the nonspatial i.i.d. unit-specific frailties in the cure model’s split and survival stages, respectively. If $V_i=W_i=0$, then Equation (4) reduces to the log-likelihood of the nonspatial “pooled” split-population survival model (without unit-specific frailties) with time-varying covariates [22,24]. Suppose, however, that unobserved unit-specific heterogeneity influences the units’ survival time or probability of entering the immune fraction (or both). Unobserved heterogeneity of this sort is addressed by incorporating the i.i.d. split-stage and survival-stage frailty terms $(V_i$ and $W_i)$ into the split-population survival model. In a Bayesian framework, the following exchangeable normal prior is employed to assess these frailties in each stage of the model:

$$W_i\sim N(0,1/\tau )\mathrm{and}{V}_i\sim N(0,1/\tau )$$

(5)

where $\tau$ is the precision parameter and each unit is specified as exchangeable to generate the prior [11,12,25].

If researchers believe that the effect of each unit-specific frailty on the unit’s risk-propensity or probability entering the immune fraction is independent of neighboring units’ frailty effects, then the nonspatial split-population survival model should be estimated with i.i.d. unit-specific frailties. It is possible, however, for the frailties to exhibit spatial heterogeneity, meaning that each unit’s propensity to be in the immune fraction as well as its survival time is influenced by unobserved factors among its neighboring units. Spatial weights can be assigned to the unit-specific frailties in the split-population survival’s split and survival stage to model this spatial autocorrelation. These spatially weighted frailties can then be incorporated via the conditionally autoregressive (CAR) approach previously developed in [26].

In the Bayesian split-population survival model, the CAR prior accounts for spatial autocorrelation in the frailties by allowing these frailities to be autocorrelated across, e.g., geographically adjacent units, where “adjacency” can be defined by the researcher based on the context. More specifically, spatial data are often represented by a lattice, in which the spatial surface is divided into a grid of units that, depending on the empirical context, can be counties, districts, countries, or other areal units. The spatially weighted frailties are then incorporated via the CAR prior by defining the relevant spatial relationship among all geographically adjacent units in an adjacency matrix $\mathbf{A}$, where each element is denoted as $a_{i{i}^{\prime }}$.

Note that $a_{i{i}^{\prime }}=1$ in $\mathbf{A}$ if units i and $i^{\prime }$ are “neighbors”. If i and $i^{\prime }$ are not neighbors, then $a_{i{i}^{\prime }}=0$. The assignment of spatial weights is incorporated into the CAR prior in order to model spatially autocorrelated frailties between adjacent units. After doing this, the frailties $V_i$ and $W_i$ are collected into the vectors $\mathbf{V}=\{V_1$,…,$V_N\}$ and $\mathbf{W}=\{W_1$,…,$W_N\}$, respectively. This facilitates the use of separate CAR priors for $\mathbf{V}$ and $\mathbf{W}$, which in turn produces the following CAR model structure:

$$\mathbf{V}|\lambda \sim \mathrm{CAR}(\lambda \left)\mathrm{and}\mathbf{W}\right|\lambda \sim \mathrm{CAR}\left(\lambda \right)$$

(6)

where $\lambda$ is the precision parameter [11,26]. The CAR$\left(\lambda \right)$ prior for $\mathbf{V}$ and $\mathbf{W}$ has a joint distribution, which is formally characterized in [14] and is described in our paper’s Supplementary Materials. The resulting conditional distributions of the spatial frailties for $\mathbf{V}$ and $\mathbf{W}$ are

$$V_i|V_{i^{\prime }\ne i}\sim N(\overline{V_i},1/\left(\lambda m_i\right)),W_i|W_{i^{\prime }\ne i}\sim N(\overline{W_i},1/\left(\lambda m_i\right)),$$

(7)

where: $\overline{W_i}=m_i^{-1}{\sum }_{i^{\prime }adji}{W}_{i^{\prime }}$; $\overline{V_i}=m_i^{-1}{\sum }_{i^{\prime }adji}{V}_{i^{\prime }}$; ${\overline{W}}_i$ and $\overline{V_i}$ are the averages of the neighboring $W_{i^{{}^{\prime }}\ne i}$ and $V_{i^{{}^{\prime }}\ne i}$, respectively, where $i^{\prime }$$adj$i denotes that $i^{\prime }$ is adjacent to i given the matrix $\mathbf{A}$; and $m_i$ is the number of these adjacencies [25,27]. Using this CAR prior approach, we can then define the log-likelihood of the spatial split-population survival model by substituting $\mathbf{V}=\left\{V_i\right\}$ and $\mathbf{W}=\left\{W_i\right\}$ in Equation (4), where ${\alpha }_{ij}=\frac{exp({\mathbf{Z}}_{ij}\gamma ,\mathbf{V})}{1+exp({\mathbf{Z}}_{ij}\gamma ,\mathbf{V})}$ is the split-stage equation.

The log-likelihood of the pooled (“nonfrailty”), nonspatial i.i.d. frailty, and spatial split-population survival models are compatible with any commonly employed parametric survival distribution. For our empirical applications, we assume a Weibull distribution for the baseline hazard, in which $\rho$ denotes the shape parameter. The density, survival function, and the hazard rate for the Weibull distribution are defined in our Supplementary Material. We use the Geweke [28] convergence test and Heidelberger and Welch [29] stationarity test in our empirical applications below to assess whether the obtained Markov chains converge to their respective stationary distributions.

2.2. Markov Chain Monte Carlo Estimation

Following standard practice for Bayesian inference [30], we assign the Multivariate Normal (MVN) prior to $\beta =\{{\beta }_1,\dots ,{\beta }_{p_1}\}$ and $\gamma =\{{\gamma }_1,\dots ,{\gamma }_{p_2}\}$, and the Gamma prior for $\rho$ with shape and scale parameters $a_{\rho }$ and $b_{\rho }$ for the Bayesian pooled (nonfrailty), nonspatial (i.i.d.) frailty, and spatial split-population parametric (Weibull) survival models:

$$\begin{array}{cc}\hfill \rho & \sim \mathrm{Gamma}(a_{\rho },b_{\rho }),\beta \sim {MVN}_{p_1}({\mu }_{\beta },{\mathsf{\Sigma }}_{\beta }),\gamma \sim {MVN}_{p_2}({\mu }_{\gamma },{\mathsf{\Sigma }}_{\gamma })\hfill \\ \hfill & {\mathsf{\Sigma }}_{\beta }\sim \mathrm{IW}(S_{\beta },{\nu }_{\beta }),{\mathsf{\Sigma }}_{\gamma }\sim \mathrm{IW}(S_{\gamma },{\nu }_{\gamma }),\hfill \end{array}$$

(8)

where $a_{\rho }$, $b_{\rho }$, $S_{\beta }$, ${\nu }_{\beta },S_{\gamma }$, and ${\nu }_{\gamma }$ are the hyperparameters in (8) and ${\mu }_{\beta }$ and ${\mu }_{\gamma }$ are random variables. Here, ${\mathsf{\Sigma }}_{\beta }$ and ${\mathsf{\Sigma }}_{\gamma }$ are estimated in a Bayesian hierarchical framework using the Inverse Wishart (IW) distribution when employing the MVN (weakly informative) prior. For Bayesian MCMC estimation of the spatial split-population parametric (e.g., Weibull) model, we additionally assign the hyperprior $p\left(\lambda \right)$ to $\lambda$ in light of the CAR prior approach. Specifically, we assign the Gamma hyperprior $\lambda \sim \mathrm{Gamma}(a_{\lambda },b_{\lambda })$ for $\lambda$ [11,12]. We specify the vague prior $(a_{\lambda },b_{\lambda })=(0.001,1/0.001)=(0.001,1000)$, as for the case of for $\rho$. To estimate the nonspatial frailty split-population survival model in this case, we assign the normal prior for the model’s split and survival-stage frailties $(V_i,W_i)$, and use the prior described above for the model’s $\beta$, $\gamma$, and $\rho$ parameters. To identify the nonspatial frailty and spatial split-population model intercepts, we impose the constraint that the frailties sum to zero, i.e., ${\sum }_i{V}_i=0$ and ${\sum }_i{W}_i=0$.

The joint posterior distribution of the Bayesian spatial split-population Weibull model with time-varying covariates is

$$\begin{array}{cc}\hfill \pi (\beta ,\gamma ,\rho ,\mathbf{W},\mathbf{V},\lambda ,{\mathsf{\Sigma }}_{\beta },{\mathsf{\Sigma }}_{\gamma }|& \mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\gamma )\propto L(\beta ,\gamma ,\rho ,\mathbf{W},\mathbf{V}|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0})\hfill \\ \hfill & \pi \left(\mathbf{W}\right|\lambda )\pi \left(\mathbf{V}\right|\lambda )\pi \left(\beta \right|{\mu }_{\beta },{\mathsf{\Sigma }}_{\beta })\pi \left(\gamma \right|{\mu }_{\gamma },{\mathsf{\Sigma }}_{\gamma })\pi \left(\rho \right)\pi \left(\lambda \right)\pi \left({\mathsf{\Sigma }}_{\beta }\right)\pi \left({\mathsf{\Sigma }}_{\gamma }\right),\hfill \end{array}$$

(9)

where the likelihood $L(\beta ,\gamma ,\rho ,\mathbf{W},\mathbf{V}|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0})$ is from Equation (4) with frailties $V_i$ collected into $\mathbf{V}=\{V_1,\dots ,V_N\}$ and $W_i$ into $\mathbf{W}=\{W_1$,…,$W_N\}$; here, $\mathbf{C}$ represents the vector of censored observations. The density, survival function, and hazard rate for this likelihood are defined in the Supplementary Materials for the Weibull case; $\pi \left(\mathbf{W}\right|\lambda )$ and $\pi \left(\mathbf{V}\right|\lambda )$ are defined via their respective conditional distributions in Equation (9), $\pi \left(\beta \right|{\mu }_{\beta },{\mathsf{\Sigma }}_{\beta })$, $\pi \left(\gamma \right|{\mu }_{\gamma },{\mathsf{\Sigma }}_{\gamma })$, $\pi \left(\rho \right),\pi \left({\mathsf{\Sigma }}_{\beta }\right)$, and $\pi \left({\mathsf{\Sigma }}_{\gamma }\right)$ are defined in Equation (8), and $\pi \left(\lambda \right)$ is the Gamma hyperprior for the spatial split-population parametric survival models. From (9), we can formally state the joint posterior distribution of the time-varying nonspatial frailty split-population parametric (Weibull) model by incorporating the frailties $V_i$ and $W_i$ defined in Equation (6) instead of $\mathbf{W}$, $\mathbf{V}$, and their respective CAR priors. The conditional posterior distributions for $\beta ,\gamma$, and $\rho$ in the pooled (nonfrailty) parametric model with time-varying covariates are

$$\begin{array}{cc}\hfill P\left(\beta \right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\gamma ,\rho )& \propto P(\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\beta ,\gamma ,\rho )\times P\left(\beta \right|{\mathsf{\Sigma }}_{\beta }),\hfill \\ \hfill P\left(\gamma \right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\beta ,\rho )& \propto P(\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\beta ,\gamma ,\rho )\times P\left(\gamma \right|{\mathsf{\Sigma }}_{\gamma }),\hfill \\ \hfill P\left(\rho \right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\beta ,\gamma )& \propto P(\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\beta ,\gamma ,\rho )\times P\left(\rho \right|a_{\rho },b_{\rho }),\hfill \end{array}$$

(10)

where $P(\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{t}\mathbf{0},\beta ,\gamma ,\rho )$ is the likelihood obtained from Equation (4) after excluding the frailty terms and $P(\beta$$|{\mathsf{\Sigma }}_{\beta })$, $P(\gamma$$|{\mathsf{\Sigma }}_{\gamma })$, and $P(\rho$$|a_{\rho },b_{\rho })$ are the priors defined in Equation (8).

The pooled (nonfrailty), nonspatial (i.i.d.) frailty, and spatial split-population survival model using the Weibull distribution can be estimated using an MCMC algorithm for Bayesian inference. To begin with, because closed form distributions for the posterior distributions of $\beta ,\gamma$, $\rho ,\lambda ,\mathbf{W}$ and $\mathbf{V}$ are not available for the spatial split-population survival model, our MCMC method’s update scheme in this case incorporates Gibbs Sampling (for estimating $\lambda$), the Metropolis–Hastings algorithm (for $\mathbf{W}$ and $\mathbf{V}$ given $\lambda$), and slice-sampling [31] for updating $\beta ,\gamma$, $\rho$. We use Gibbs sampling for $\lambda$, as it is easier to sample from the conditional distribution (which is known) in this case and because the joint distribution is not known explicitly. We employ the Metropolis–Hastings algorithm for $\mathbf{W}$ and $\mathbf{V}$ given $\lambda$ because it is difficult to sample from the conditional distribution for $\mathbf{W}$ and $\mathbf{V}$. Finally, we use slice-sampling for updating $\beta ,\gamma$ and $\rho$, as it requires little tuning because the slice width adapts quickly to the distribution and sampler performance. Furthermore, considering that slice-sampling draws from the posterior samples from any prior distribution as long as these distributions have a reasonable value range of parameters, this sampling algorithm provides researchers with flexibility in the choice of prior distribution. This permits the use of informative, weakly informative, or non-informative priors.

The MCMC algorithm described above proceeds as follows:

• Choose a starting point ${\beta }_0,{\gamma }_0$, ${\rho }_0$, ${\lambda }_0$ and corresponding ${\mathbf{W}}_0=\{W_1,\dots ,W_N\}$ and ${\mathbf{V}}_0=\{V_1,\dots ,V_N\}$, then set $k=0$.
• Update ${\mathsf{\Sigma }}_{\beta }\sim \pi \left({\mathsf{\Sigma }}_{\beta }\right|\beta )$, ${\mathsf{\Sigma }}_{\gamma }\sim \pi \left({\mathsf{\Sigma }}_{\gamma }\right|\gamma )$, $\lambda \sim \pi \left(\lambda \right|\mathbf{W},\mathbf{V})$ using Gibbs sampling. The closed form of the full conditional distributions for $\pi \left({\mathsf{\Sigma }}_{\beta }\right|\beta )$, $\pi \left({\mathsf{\Sigma }}_{\gamma }\right|\gamma )$, $\pi \left(\lambda \right|\mathbf{W},\mathbf{V})$ are derived and defined in the Supplementary Materials.
• Update $\beta \sim \pi \left(\beta \right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{W},\mathbf{V},\gamma ,\rho ,{{}^{–}}_{\beta },{\mathsf{\Sigma }}_{\beta })$, $\gamma \sim \pi \left(\gamma \right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{W},\mathbf{V},\beta ,\rho ,{{}^{–}}_{\gamma },{\mathsf{\Sigma }}_{\gamma })$, and $\rho \sim \pi \left(\rho \right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{W},\mathbf{V},\beta ,\gamma ,a_{\rho },b_{\rho })$ using the slice sampler with stepout and shrinkage (Neal, 2003); see the Supplementary Materials for details on performing the slice sampling operation in this step.
• Update $\mathbf{W}\sim \pi \left(\mathbf{W}\right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{V},\beta ,\gamma ,\rho ,\lambda )$ and $\mathbf{V}\sim \pi \left(\mathbf{V}\right|\mathbf{C},\mathbf{X},\mathbf{Z},\mathbf{t},\mathbf{W},\beta ,\gamma ,\rho ,\lambda )$ via Metropolis–Hastings.
• Set $k=k+1$, then return to Step 2 and repeat for K iterations.

The MCMC algorithm for estimation of the nonspatial frailty model is similar to the steps delineated above, with the exception of the nonspatial i.i.d. frailties $V_i$ and $W_i$ in this model being updated via Metropolis–Hastings with the proposal variance as the conditional prior variance for these frailties. To estimate the pooled (nonfrailty) model, we use the following MCMC algorithm, as closed forms for the posterior distributions of $\rho$, $\beta$, and $\gamma$ are not available:

• Choose the initial values of $\beta ,\gamma$, and $\rho$, then set $m=0$.
• Update ${\mathsf{\Sigma }}_{\beta }$ and ${\mathsf{\Sigma }}_{\gamma }$ via Metropolis–Hastings; see the Supplementary Material for the closed form of the full conditional distributions for ${\mathsf{\Sigma }}_{\beta }$ and ${\mathsf{\Sigma }}_{\gamma }$.
• Update $\beta$, $\gamma$, and $\rho$ using the slice sampler with stepout and shrinkage, as described in the Supplementary Materials.
• Repeat Steps 2 and 3 until the chain converges.
• After M iterations, summarize the parameter estimates using posterior samples.

3. Monte Carlo Simulations

We conducted three Monte Carlo (MC) experiments to compare the performance of the nonspatial SP Weibull models with and without i.i.d. frailties to our spatial split-population Weibull model. The design of our MC experiments and the results from these experiments are presented in greater detail in the Supplementary Materials. We focus on the Weibull case here, as the empirical applications below use the Weibull survival distribution; however, our Monte Carlo simulation results hold for other parametric (e.g., log-logistic) distributions as well.

More specifically, our MC experiments simulate a split-population Weibull distributed outcome variable that exhibits spatial autocorrelation across neighboring units in each stage. For all experiments, we consider sample sizes of $N=100$, 400, 1000, 1500, and 2000. Note that $N=100$ and $N=400$ respectively correspond to a small and moderate sample size, while $N=1000,1500$, and 2000 represent a relatively larger sample. For each model in our MC experiments, we include one survival-stage covariate ${\mathbf{x}}_1$ and two split-stage covariates ${\mathbf{z}}_1$ and ${\mathbf{z}}_2={\mathbf{x}}_1$, as the same covariate may be included in both stages. We incorporate information about the spatial relationship between units in our simulated data via an adjacency matrix $\mathbf{A}$. To generate $\mathbf{A}$, we consider a hypothetical space with five areal units (e.g., countries), with each unit having at least one adjacent “neighbor.” This spatial relational information is then incorporated into the simulated data generation process, which follows an SP Weibull distribution (see Supplementary Materials for details).

Next, recalling that the split-stage equation in the spatial split-population survival model is provided by a binary response function that captures the effect of covariates ${\mathbf{Z}}_i$ and the associated parameter vector $\gamma$ on the probability of units entering the immune fraction ($\alpha$), we have a case in which the more likely a greater share of units is to enter the immune fraction, the higher the immune fraction level. Hence, for the MC experiments, we set the true $\gamma$ values that affect the immune fraction (via $\alpha$)—calculated as the mean value of the binary response function ${\alpha }_i=\frac{exp({\mathbf{Z}}_i\gamma +V_i)}{1+exp({\mathbf{Z}}_i\gamma +V_i)}$) for all i in our N-sample data—using the pre-set true $\gamma$ value and the randomly generated variables ${\mathbf{Z}}_i$ (as well as the $\mathbf{V}$ spatial frailties for the spatial d.g.p. This permits us to adjust the immune fraction level in a way that is consistent with the model’s split-stage. Finally, for each experiment, we use 500 iterations in the MCMC, 100 burn-ins, and a thinning of 1, and assess the convergence of the Markov chain via trace-plots and the Geweke convergence test.

Using these experimental conditions, we now turn to assessing our nonspatial and spatial models of interest for three experiments. In the first MC experiment, we compare the performance of our nonspatial and spatial split-population Weibull models when the fraction of the immune subpopulation is fixed at 25% and the proportion of units that share spatial frailties is held at 40%. The results from this MC experiment reveal that our spatial split-population Weibull model outperforms both the nonspatial split-population Weibull models in retrieving the true theoretical values of the split-stage ($\gamma$) and survival-stage ($\beta$) covariates (Figure 1) along with the spatial frailties in both stages (Figure S1, Supplementary Material) for small and moderate sample sizes as well as for the relatively larger sample sizes listed above. Thus, the spatial split-population Weibull model should be favored over the nonspatial split-population models even when a low number of observations have spatially-dependent frailties in a split-population survival framework. The trace plots show stability, and the models pass the Geweke convergence test.

In the second MC experiment, we compare the nonspatial and spatial split-population Weibull models’ performance when the immune fraction of the simulated split-population Weibull-distributed outcome variable remains at 25% and the proportion of units that share spatial frailties is 30%, 40%, 60%, and 80%. Results from Figure 2, which displays the mean RMSEs for $\widehat{\beta }$ and $\widehat{\gamma }$ for $N=100,400,1000,1500,$ and 2000, demonstrate that the spatial split-population Weibull model substantially outperforms the nonspatial split-population Weibull models at all levels of spatial autocorrelation; the mean RMSEs for the spatial SP Weibull models are always negligible (close to 0) for all the sample sizes that are examined in the experiment, while those of the two nonspatial SP Weibull models are critically high, indicating a considerable level of bias. Moreover, the results show that the nonspatial split-population Weibull models’ split and survival-stage covariates and nonspatial frailties exhibit deteriorating coverage and lower efficiency at all levels of spatial autocorrelation. By contrast, the spatial split-population Weibull model recovers the true theoretical values of the split and survival-stage covariates and the spatial frailties in both stages of the model with more accuracy, coverage, and efficiency. Therefore, for the various sample sizes considered here, the spatial split-population Weibull model outperforms the nonspatial models at all levels of spatial autocorrelation. Again, the trace plots show stability, and the models pass the Geweke convergence test.

The third MC experiment applies the d.g.p. from Experiment 1, where the share of units with spatially dependent frailties is held at 40% and the size of the immune fraction is varied from 25% to 33%, 40%, 48%, and finally 60%. Figure 3a,b, which additionally shows the mean RMSEs for $\widehat{\beta }$ and $\widehat{\gamma }$ with $N=100,400,1000,1500,$ and $2000$, reveals that the split and survival stage results strongly favor the spatial split-population Weibull model over the other two nonspatial split-population Weibull models for all the sample sizes examine here and at all immune fraction levels; the mean RMSEs are close to 0, and are negligible in the spatial SP Weibull models, while those for the two nonspatial SP Weibull models are considerably higher. This indicates that the $\widehat{\beta }$ and $\widehat{\gamma }$ parameters are biased in the latter models. Furthermore, the retrieved values of the spatial frailties from the spatial split-population Weibull model rapidly converge to their true values with high coverage probabilities. Thus, if the true d.g.p. is split-population Weibull in which the unit-specific frailties exhibit spatial autocorrelation, then the spatial split-population Weibull should be estimated when the size of the immune fraction is 25% or above for small, moderate, and even relatively large sample sizes.

We conducted three additional MC experiments that we do not report here in order to save space; these are briefly presented in the Supplementary Materials (see Figure S2 and Tables S3 and S4). Briefly, in one set of MC experiments (Experiment 4) we increased both the immune fraction and the share of units with spatially dependent frailties. In another experiment (Experiment 5), we re-evaluated our primary MC results using an alternative prior. In the third set of MC experiments (Experiment 6), we compared our model’s performance to a Bayesian spatial SP model that incorporates spatial frailties in just the survival stage. For each of these additional experiments, we set $sims=100$ and evaluated model performance for a variety of different sample sizes. These additional MC experiments revealed that, unlike the nonspatial models, both the retrieved values of the split and survival-stage parameters and the spatial frailties from the spatial SP Weibull model converge to their true values with high coverage probabilities when the true d.g.p. is SP Weibull in which the unit-specific frailties exhibit spatial autocorrelation.

4. Empirical Applications

Our MC experiments suggest that if spatial autocorrelation exists in the true d.g.p., then failing to account for it leads to faulty inferences. In applied settings, this means that if there are a priori theoretical reasons to suspect that the survival times and immune fractions of interest are spatially clustered, where our spatial frailty approach is superior to nonfrailty or i.i.d. frailty split-population models. Below, we apply our Bayesian spatial split-population Weibull model to survival data from two published studies in Political Science about (1) democratic regime survival [17] and (2) the duration of post-civil war peace (i.e., before civil war recurs) [18]. In both cases, we discuss theoretical reasons and empirical evidence suggesting that these processes exhibit spatial clustering, then compare results from our Bayesian spatial frailty model to those of nonspatial specifications.

4.1. Democratic Consolidation and Survival

Comparative political scientists often conceptually distinguish between transitional democracies, which can revert to authoritarian rule, and consolidated democracies, in which the “democratic regime becomes sufficiently durable that democratic breakdown is no longer likely” [32] (p. 743). Previous empirical analyses of democratic consolidation, which typically employ discrete choice models, consistently find that wealth (measured by GDP per capita) has a positive and highly significant effect on the probability of democratic consolidation [32,33]. On the other hand, presidential systems (as opposed to parliamentary systems) have a negative but weakly significant or insignificant impact on democratic consolidation, while the association between past authoritarian institutions and democratic consolidation is inconsistent [32,34]. Research on the survival of democratic regimes (which often uses conventional parametric survival models) usually finds that, unlike presidential systems, economic growth and parliamentary systems help democracies to endure [35,36], while democracies preceded by military rule revert to dictatorships more quickly [34,37]. The most consistent finding, however, is that GDP per capita has a strong positive influence on survival of democracy, leading a number of scholars to infer that democratic survival “increases monotonically with per capita income” and then endures indefinitely after GDP per capita reaches approximately USD 6000 [38] (p. 165).

Although these insights are important, Ref. [17] emphasizes that by employing standard duration or discrete choice models these studies assume that all democracies face the same baseline risk of reversal to authoritarianism. This assumption is unjustified, as the population of democracies includes “at-risk” transitional democracies along with an “immune fraction” of fully consolidated democracies for which the risk of authoritarian reversal is negligible. Hence, the observed survival of democracy results from two separate processes: “democracies that survive because they are consolidated and those democracies that are not consolidated but survive because of some favorable circumstances” [17] (p. 153). Considering these two subpopulations, Ref. [17] then re-examines extant findings about democratic consolidation and survival by estimating via MLE parametric (Weibull) nonspatial split-population survival models (with and without i.i.d. frailties) on a dataset of democratic spells across 133 countries between 1789 and 2001. All right-censored observations in his data are either consolidated or transitional democracies that have not yet reverted to authoritarian rule. Thus, the split-stage in his SP survival model estimates the probability of democratic consolidation (62% of his cases are right-censored), while the survival stage estimates the duration of democracy among cases that eventually experience an authoritarian reversal. He incorporates seven covariates in both stages: GDP per capita, GDP growth, Presidential and Parliamentary systems, and previous Military, Civilian, and Monarchical dictatorships.

Briefly, [17] found that GDP per capita has a positive and statistically significant effect on the probability of democratic consolidation in the split-stage, while presidential systems and democracies preceded by military dictatorships are less likely to consolidate. Economic growth helps “transitional” democracies to survive longer, though there is no statistically significant relationship between GDP per capita and democratic survival among these at-risk regimes. Most other covariates are statistically insignificant, though the insignificant coefficient for presidential systems is notably positive in the survival stage.

Despite this important contribution to the extensive literature on the durability of democracies, the nonspatial split-population survival models employed in the original study assumed that neither the likelihood of democratic consolidation nor the prospects for democratic survival exhibit spatial autocorrelation. This assumption may be untenable, as democracies tend to cluster in space. Indeed, “since 1815, the probability that a randomly chosen country will be a democracy is about 0.75 if the majority of its neighbors are democracies, but only 0.14 if the majority of its neighbors are nondemocracies” [39] (p. 916). Research in political science has found that geographical proximity to (consolidated) democracies not only encourages democratic transition in authoritarian regimes, it increases the odds of consolidation and survival of nascent democracies, as stable democracies create a “regional production chain” of democratic institutions, practices, and norms that are conducive to democracy [40] (p. 25; and see [32,39,41]). Thus, democratic clustering reinforces democratic norms, making it costly for elites to engage in democratic backsliding [41]. In this way, democratic neighborhood effects may have important latent influences on both democratic reversal and democratic survival.

Although there are clear theoretical reasons to expect spatial autocorrelation in the d.g.p. of democratic regime survival and consolidation, we can additionally use common tools in spatial statistics to diagnose the extent of spatial clustering in each year of the data. Specifically, we conduct two pre-estimation tests by calculating (1) the join count and (2) the Global Moran’s I statistic for each cross-section of democracies in the data. The join count is a measure of the extent to which the number of observed areal units that are adjacent and of the same category is greater than or less than what is expected if the spatial distribution of those categories were random [42]. In general, in a setting with two discrete categories A and B, the join count test statistic is

$$Z\left(AB\right)=\frac{AB-E\left(AB\right)}{\sqrt{{\sigma }_{AB}^2}},$$

(11)

where AB and E$\left(AB\right)$ are the observed and expected counts of adjacent units in categories A and B, respectively, and ${\sigma }_{AB}^2=E\left(A{B}^2\right)-E{\left(AB\right)}^2$. Positive statistics indicate spatial dispersion (units of the same category are further from each other than expected by chance), while negative statistics indicate positive spatial clustering (units of the same category are more likely to be adjacent than what is expected by chance.

For this application, we construct a separate cross-sectional adjacency matrix with elements $a_{i{i}^{\prime }}$ for each year in the data, wherein proximate pairs of democratic countries (within 800 km of each other) are assigned a weight of 1 ($a_{i{i}^{\prime }}=1$). Our outcome of interest for the join count analysis is whether a country is identified as being “at risk” of democratic reversal in the original dataset.

In addition, we use the Global Moran’s I statistic to assess the number of years that democratic regimes survive as a group of clusters in space. Global Moran’s I is an inferential statistic that measures the direction and degree of spatial clustering in continuous data [43]. Positive statistics indicate positive spatial clustering of similar values of the continuous variable of interest, while negative statistics indicate that dissimilar values are more likely to be proximate than if they were distributed randomly in space. Using the same adjacency matrix described above, we use the Global Moran’s I to evaluate whether democracies that have survived for similar periods of time exhibit spatial heterogeneity (clustering or dispersion).

We report the join count and Moran’s I tests in detail in the Appendix A (Figure A1) provided at the end of the paper. Briefly, the results clearly indicate significant spatial clustering in both the probability of democratic reversal and the survival rates of democracies, particularly in the post-World War 2 period. Next, we replicate the above analysis using our Bayesian spatial SP survival model in order to compare our results to the original nonspatial models with and without i.i.d. frailties. Because the original analysis of these data used maximum likelihood estimation, we used the same for the nonspatial frailty and nonfrailty models in order to exactly replicate the previous results. Our spatial SP Weibull model incorporates spatially-weighted frailties across neighboring democracies via the adjacency matrix A. We construct a matrix A with elements $a_{i{i}^{\prime }}$ such that $a_{i{i}^{\prime }}=1$ for each year if the capital of country i is less than 800 km from the capital of country $i^{\prime }$ and $a_{i{i}^{\prime }}=0$ if countries i and $i^{\prime }$ are greater than 800 km from each other. Using geographic proximity as the spatial relationship of interest is appropriate, as it allows the frailties to be correlated with those of neighboring democracies rather than assuming spatial independence even within the same regions. Considering our Bayesian MCMC estimation approach, we incorporate the spatial information in A by employing separate CAR priors for the frailty terms vector V (split-stage) and W (survival-stage), which implies a CAR structure of $\mathbf{V}|\lambda$∼ CAR$\left(\lambda \right)$ and $\mathbf{W}|\lambda$∼ CAR$\left(\lambda \right)$. The spatial SP Weibull model is estimated based on the sample from [17] using the MVN prior and our MCMC algorithm described earlier and assigning the Gamma hyperprior for $\lambda$. Here, we use the hyperparameters $a=1$, $b=1$, $S_{\beta }$ = $I_{p1}$, $S_{\gamma }$ = $I_{p2}$, ${\nu }_{\beta }=p1$, ${\nu }_{\gamma }=p2$. Recall that ${\mathsf{\Sigma }}_{\beta }$ is the variance of the MVN prior of the vector $\beta$ for p1-dimensional survival stage covariates and that ${\mathsf{\Sigma }}_{\gamma }$ is the MVN’s prior of the vector $\gamma$ for p2-dimensional split-stage covariates. Hence, when we employ the Inverse Wishart (IW) distribution to estimate both ${\mathsf{\Sigma }}_{\beta }$, in which ${\nu }_{\beta }$ is the hyperparameter, and ${\mathsf{\Sigma }}_{\gamma }$, in which ${\nu }_{\gamma }$ is the hyperparameter, we adopt the values p1 for ${\nu }_{\beta }$ and p2 for ${\nu }_{\gamma }$. Finally, $\lambda \sim \mathrm{Gamma}(a_{\lambda },b_{\lambda })$ with a vague prior $(a_{\lambda },b_{\lambda })=(0.001,1/0.001)$. Our Bayesian SP Weibull model results are based on a set of 50,000 iterations after 4000 burn-in scans and thinning of 10.

We begin our analysis of the split-stage results by examining choropleth maps (Figure 4a,b) that illustrate the posterior means of the spatial frailties obtained from the spatial SP Weibull model. The split-stage map (Figure 4a) reveals that there are distinct spatial bands in the frailties, which range from $-0.725$ to $0.716$ with a corresponding standard deviation of 0.31. The spatial patterns in the map suggest that there is strong spatial clustering in the underlying factors linked to democratic consolidation, as states with a higher baseline risk for democratic consolidation are in similar geographic neighborhoods, whereas those with lower propensities are clustered in separate regions.

Figure 5 displays the results for each covariate from the replicated models and our spatial frailty model. For the nonspatial models, the points represent coefficient estimates and the bars represent 90% confidence intervals. For admittedly rough comparability purposes, the dots in the figure represent posterior means for the spatial frailty models, while the bars represent symmetric 90% credible intervals. Although these are certainly not perfect comparisons, our goal here is simply to illustrate the applicability of our model and the way in which accounting for spatial autocorrelation can affect inferences.

The dot-whisker plots in Figure 5 show that while Svolik’s results for GDP growth, Monarchy, and Civilian are similar across the nonfrailty and nonspatial frailty SP Weibull models, the differences in the split-stage results between these nonspatial models and our spatial frailty model are more pronounced. For instance, the Presidential and Military covariates are each negative and highly significant in the nonspatial SP Weibull models with and without i.i.d. frailties. By contrast, the negative estimates for Presidential and Military are each highly unreliable in the spatial SP Weibull model’s split-stage equation. Thus, after we explicitly account for spatial autocorrelation in the split-stage of the SP survival model, the relationship between each of the covariates noted above and the probability of democratic consolidation is considerably attenuated. In nonspatial SP Weibull models from [17], the Parliamentary dummy’s split-stage estimate is positive, albeit insignificant. However, the estimate of Parliamentary in the spatial SP Weibull model’s split-stage is negative (though not reliably so). Thus, the results from our spatial model raises doubts about prior claims that parliamentary systems are strongly associated with democratic consolidation.

Finally, we consider the split-stage parameter estimate for GDP per capita, which is positive and statistically significant in the original nonspatial SP Weibull models with and without i.i.d. frailties. In contrast, the split-stage estimate of GDP per capita is negative (though insignificant) in our spatial SP Weibull model. Hence, the widely accepted positive association between higher per capita income and democratic consolidation is neither consistent nor robust when accounting for spatial autocorrelation among neighboring democracies.

Turning to the survival stage results, we first consider the choropleth map in Figure 4b, which illustrates posterior means of the spatial frailties obtained from the Spatial SP Weibull model. The spatial frailty values vary from $-0.95$ to $0.859$, with a corresponding standard deviation of 0.313. These maps again reveal spatial clustering associated with democratic regime survival; those democracies with greater underlying propensity for democratic survival are located near countries with similar propensities, while those with a lower propensity for democratic survival are located in disparate geographic areas.

The plots in Figure 6 reveal additional differences between the original nonspatial split-population model results and the new results from the Bayesian spatial split-population model. For instance, although the original study found that the survival stage estimate of Monarchy is positive and highly significant in the nonspatial models, this relationship is insignificant and negative in the spatial SP Weibull model. Hence, the association between democracies that were previously ruled by a monarch and democratic durability is tenuous after accounting for the influence of spatial autocorrelation on democratic survival. Next, the original survival stage estimate for Military in both the nonspatial SP Weibull models is negative (though insignificant). In the spatial SP Weibull model, Military is again negativ; however, unlike the nonspatial models, in this case it is statistically reliable. This suggests that not accounting for neighborhood democracies can lead researchers to underestimate the relationship between democratic durability and democratic states that were preceded by military rule. Finally, while [17] challenged the confidence of previous findings with respect to the relationship between GDP per capita and democratic survival (e.g., [35,36]), the influence of per capita income on democratic survival in our spatial SP Weibull model is positive and statistically reliable, consistent with the previous literature. Taken together, a re-examination of these data on democratic survival and consolidation using the spatial SP survival model leads to new inferences about important political phenomena.

4.2. Post-Civil War Peace Duration

To further demonstrate the applicability of our framework, in this section we use our spatial SP survival model to re-investigate previous findings that suggest information transparency and other political freedoms can increase post-civil war peace survival [18]. The normative importance of consolidating peace after civil war has motivated a wide body of research in Political Science on why certain civil wars recur and others do not. Much of this literature has focused on the characteristics of the initial war and termination [44,45] and the characteristics of the post-war environment, including the presence of third-party or U.N. intervention [46,47,48]; ref. [18] contributes to this literature by arguing that increased political accountability and civil liberties in the post-war period can augment the costs of reneging on an agreement and allow governments to credibly commit to not resuming violence. Her primary testable expectation is that “civil wars that are fought against governments with limited accountability should be more likely to repeat themselves than civil wars in countries with highly accountable governments” [18] (p. 1248). Estimating conventional survival models using a sample of 77 post-civil war peace spells during 1945–2009, Ref. [18] finds support for her hypotheses in various measures of civil liberties, democratic institutions, and rule of law; however, one of her primary variables of interest, Press Freedom, does not reach standard frequentist levels of statistical significance.

The conventional survival models used in the aforementioned study, however, assume that all wars recur at some point. A split-population survival model may be more appropriate for studying post-war peace survival, as certain wars are at risk for recurring while others are structurally distinct cases in which peace is “consolidated” or one side is eliminated entirely, in which the same conflict cannot recur. Moreover, there are theoretical and empirical reasons to believe that spatial autocorrelation may influence both the risk of peace consolidation or of renewed war in these data. Previous research has shown that diffusion processes may lead to conflict contagion (e.g., [13,49,50]), while peace stability is regionally clustered due to the clustering of observable and unobservable political attributes (e.g., [39,51,52]). In the split stage, post-war countries might never return to violence if they are surrounded by similarly peaceful countries that have demonstrated the institutional capability and political interest to prevent the resurgence of violence in the region [49], or if latent localized interests among elites or civilian populations to contain violence—for instance, to prevent war recurrence in their own contiguous countries—are clustered in space [53]. In the survival stage, stable institutions [51] or other latent geographic factors that transcend the borders of a single state may influence the time that rebels take to remobilize or that governments take to re-engage with dormant dissident movements. If any of these unobserved factors are not randomly distributed in space, then failing to account for this heterogeneity will lead to faulty inferences [13].

As in the previous application, we conducted Moran’s I and join count tests to evaluate whether spatial autocorrelation exists between peace survival rates and countries that experience a civil war and never experience another. The results (see Figure A2 in the Appendix A) clearly indicate spatial clustering and dispersion, especially after 1960 and prior to 1900. This pre-estimation empirical evidence suggests that spatial autocorrelation should be taken into account when modeling the survival and consolidation of post-civil war peace.

In our analysis of post-civil war peace, we use replication data from the aforementioned study and specify an SP model with four covariates in the split-stage: Press Freedom; whether the previous conflict ended in an outright Victory; percentage of Mountainous terrain; and GDP/capita. In the survival stage, we include our main variable of interest, Press Freedom, and control for GDP/capita, whether the previous conflict ended in a Peace Agreement, the Intensity of the previous conflict, Ethnic Factionalization, the presence of UN Peacekeeping forces, whether the previous conflict was over Territory, whether the country has some Non-Contiguous territory that can be used to expedite rebel mobilization, and the percentage of Mountainous terrain. We focus on Press Freedom because it closely captures information transparency, which is an important mechanism of interest in the original theory, and because the results for this variable were inconsistent with the other findings in the study.

To examine the effects of these covariates on both the probability of peace consolidation and the survival of peace, we estimate a nonfrailty SP Weibull model and a spatial SP Weibull model. Following the original analysis of these data, we estimate the nonfrailty model using maximum likelihood. For the spatial model, we again define spatial proximity as $a_{i{i}^{\prime }}=1$ if the distance between the capitals of state i and $i^{\prime }$ is less than 800 km and $a_{i{i}^{\prime }}=0$ otherwise [54], though we found similar results when increasing these distances to 2000 and 2500 km. We allow the frailties between neighboring units to be spatially correlated by employing separate CAR priors for the frailty vectors $\mathbf{V}$ and $\mathbf{W}$, which implies a CAR structure of $\mathbf{V}|\lambda$∼ CAR$\left(\lambda \right)$ and $\mathbf{W}|\lambda$∼ CAR$\left(\lambda \right)$. The spatial SP Weibull model is estimated using the multivariate normal prior. For the slice-sampling (MCMC) algorithm, we specify the hyperparameters as $a=1$, $b=1$, $S_{\beta }$ = $I_{p1}$, $S_{\gamma }$ = $I_{p2}$, ${\nu }_{\beta }=p1$, and ${\nu }_{\gamma }=p2$, and assign the Gamma hyperprior $\lambda \sim \mathrm{Gamma}(a_{\lambda },b_{\lambda })$ for $\lambda$ with vague prior $(a_{\lambda },b_{\lambda })=(0.001,1/0.001)$. We estimated the model with 50,000 iterations and 38,000 burn-ins. Every parameter passes the Geweke [28] convergence test and Heidelberger and Welch [29] stationarity test.

Turning first to the choropleth maps of the spatial frailty values from the spatial SP Weibull model in Figure 7a,b, the split-stage frailties (V) range from −1.19 to 1.49 with a standard deviation of 0.5503, and the survival-stage frailties (W) range from −1.21 to 0.727 with a standard deviation of 0.4498. In both stages, there seem to be regional clusters of frailty values.

The $\beta$ and $\gamma$ results for both the nonspatial and spatial SP Weibull models are reported in Figure 8 and Figure 9. The dots and bars are interpreted in the same way as in the previous application. Much like the original analysis, the nonspatial SP model reveals no clear evidence that information transparency affects the survival of peace after conflict; the coefficient estimates for Press Freedom are statistically insignificant in both stages. In fact, the only (weakly) significant result in the split-stage of the nonspatial model is GDP/capita, indicating limited evidence that increased economic prosperity can increase the probability of a civil conflict never recurring. Mountains and previous Victory appear to have no relationship with peace consolidation, contrary to extant findings (e.g., [44]). The survival stage results of the nonspatial Weibull model reveal that conflicts over Territory and Ethnic Factionalization are associated with longer post-conflict peace periods, while Non-Contiguous territory and previous war Intensity are associated with shorter peace periods, at least among countries at risk of recurrent conflict.

When we account for spatially autocorrelated frailties, however, we find remarkably different results. First, although conflicts over Territory remain positively associated with post-war peace survival, we cannot reliably conclude that Non-Contiguous territory, Intensity, or Ethnic Factionalization have systematic relationships with peace in at-risk conflict locations. In addition, we find that better economic opportunities are no longer reliably associated with post-war peace consolidation, though they do have a reliable statistical relationship with the survival of peace in at-risk countries. This finding is important because it suggests that while a stronger economy can make it more difficult for rebels to re-mobilize, it is not a panacea for permanently exiting the “conflict trap,” as previous research has suggested [55].

Finally, although Press Freedom appears to have no real relationship with post-war peace survival among at-risk cases, it is positively and reliably associated with the consolidation of peace after war. In other words, although a country having a higher level of press freedom may not elongate a temporary peace when conflict is likely to recur, press freedom is associated with a significant decrease in a country’s overall susceptibility to renewed war. Although there remains considerable debate over whether democratization after a conflict hinders or helps peace to endure (e.g., [56]), this evidence suggests that information transparency and an independent media may lead to a far more durable, if not permanent, period of domestic stability after a conflict ends.

5. Discussion and Conclusions

This has article developed a parametric spatial Split-Population Survival model that accounts for both the probability that some observations are immune to an event of interest and the tendency of underlying risk factors associated with political processes to cluster in space. While the model builds upon previous work [11,12,57], it is unique in that it allows for spatially autocorrelated frailties in the split and survival stages. The model incorporates time-varying covariates in both the split and survival stages. These features allow researchers to explicitly model and statistically account for spatial heterogeneity that may influence the probability of an observation becoming immune from an event as well as the duration of a process among units considered “at-risk” in panel survival data. Our innovation is distinct from extant works on spatial statistics that address different types of spatial dependence in settings with continuous or binary dependent variables or in conventional survival models [12,58,59], as it relaxes the assumption that all observations eventually experience the event of interest. Our MC experiments reveal that, unlike nonspatial models, our spatial split-population survival model provides accurate estimates when SP survival data exhibit spatial autocorrelation. Finally, we apply our model to previously published data on widely studied phenomena in political science in the contexts of democratic regime survival and the durability of post-civil war peace. After accounting for the immune fraction and spatial autocorrelation in these applications, we find evidence contrary to previous studies’ original findings, particularly in the first-stage “cured” fraction equation.

Future work can build upon our model to devise estimation routines for survival data with multinomial outcomes (e.g., competing risks), recurrent events, or variable time trends; a particularly useful next step would be to apply this approach to a semi-parametric setting [3,60]. Extant work on nonparametric solutions could benefit from this approach [61]. Although our contribution is significant in Bayesian survival modeling by allowing for spatial frailties in the split stage rather than only in the survival stage (e.g., [11]), future extension could increase the flexibility of our model by incorporating spatial frailties only in the estimation of the immune fraction, while the survival stage includes only i.i.d. frailties or no frailties. This could be useful in applied cases where a researcher believes that the processes leading to immunity are spatially clustered while the survival probabilities of those at risk of the event are not.

Of course, any Bayesian analysis can be sensitive to the choice of the prior distribution; thus, applied researchers should take care to not limit themselves to the priors implemented here depending on their topic of study. Future analytical work could continue to investigate our model’s performance using alternative priors, such as the Gamma distribution [62]. Furthermore, further application of the model could benefit from the development of goodness-of-fit tests for distributional assumptions in the data [63] and additional tools for model specification [60].

Beyond the applications presented here, this Bayesian spatial split-population survival model can be a useful tool for researchers interested in anything from success rates of vaccines or other medical treatments, to customer analytics around the purchase of new products, to coups d’état around the world. All of these phenomena and more tend to involve some fraction of the population under investigation being probabilistically immune from the event for some measurable reason(s) while being influenced by the spatial clustering of unobserved yet meaningful factors.