Power laws in the Roman Empire: a survival analysis.

The Roman Empire shaped western civilization, and many Roman principles are embodied in modern institutions. Although its political institutions proved both resilient and adaptable, allowing it to incorporate diverse populations, the Empire suffered from many conflicts. Indeed, most emperors died violently, from assassination, suicide or in battle. These conflicts produced patterns in the length of time that can be identified by statistical analysis. In this paper, we study the underlying patterns associated with the reign of the Roman emperors by using statistical tools of survival data analysis. We consider all the 175 Roman emperors and propose a new power-law model with change points to predict the time-to-violent-death of the Roman emperors. This model encompasses data in the presence of censoring and long-term survivors, providing more accurate predictions than previous models. Our results show that power-law distributions can also occur in survival data, as verified in other data types from natural and artificial systems, reinforcing the ubiquity of power-law distributions. The generality of our approach paves the way to further related investigations not only in other ancient civilizations but also in applications in engineering and medicine.

Introduction

As complementation to continuing studies by historians, anthropologists and social scientists, ancient civilizations can be analysed quantitatively, considering mathematical and computational modelling [1,2]. Ancient civilizations are examples of a complex system composed of agents that interact, collaborate, and compete for power and resources. One example of such a system is the Roman Empire, which started with Augustus (d. 14 CE) and ended with Romulus Augustulus, with the Germanic invasion from the North in 476. After this, the eastern Roman Empire, with Constantinople as its capital, continued to exist until 1453 [3]. The Roman Empire influenced western civilization and contributed to many aspects of our culture [4]. Although the first two centuries of the Empire were characterized by stability and prosperity, known as the Pax Romana, the Roman Empire underwent several crises, including many violent deaths of Roman emperors. A recent study suggested that among the 69 rulers, 43 emperors suffered a violent death [5], including homicide, suicide or death in combats with a foreign enemy of Rome. These conflicts produced patterns in the length of time that can be identified by statistical analysis, as we show in this paper.

We propose a power-law model to describe the time-to-violent-death, which measures the time before a Roman emperor suffers a violent death since the beginning of his reign. This variable is interesting from the historical point of view as it can suggest some patterns that governed the dynamics of the Roman Empire. For example, if deaths occur independently, without any defined rule, we should observe an exponential distribution of this survival time, where the hazard function is constant. Otherwise, some other distribution should be more suitable to model the time-to-violent-death. Recently, Saleh [5] analysed the time-to-violent-death of the 69 emperors of the unified Roman Empire, from Augustus (d. 14 CE) to Theodosius (d. 395 CE). The author considered statistical tools of survival data analysis to model their time-to-violent-death and verified that a mixture of two Weibull distributions captures this variable well. However, as we show here, that model does not seem to be easily extended to the long-term perspective observed for some Roman emperors who did not face a violent death. Moreover, the dataset analysed covers only part of the Roman emperors, leaving out, for instance, the Byzantine Empire. To address these limitations, we show that a power-law distribution with change points can better model the time-to-violent-deaths. This power-law behaviour in time distribution has been observed in many complex systems, including bacterial persistence [6], earthquakes [7-9], in solar radiophysics [10], stock price fluctuations [11] and tree-limb branching [12]. The cause of this distribution is related to the presence of extreme events, such as large earthquakes. In the case of the Roman Empire, we observe that an extreme event is an emperor surviving for many years since the beginning of his reign. Our model shows that most of Rome’s emperors suffered violent deaths during the initial years of their reign.

We study the distribution of time-to-violent-death based on survival analysis, which refers to a set of statistical approaches used to investigate the time it takes for an event of interest to occur [13]. The term survival analysis originates from clinical research, where predicting the time to death, i.e. survival, is often the primary goal. The maximum-likelihood estimator (MLE) is used to estimate the best distribution for all the considered emperors, considering data in the presence of censoring and long-term survivors. Indeed, the inclusion of long-term survival enables a more accurate prediction regarding the emperors who did not have a violent death than in the study by Saleh [5]. By considering the dataset by Saleh [5], we also verify that the risk after assuming the throne is very high for a Roman emperor, but this risk systematically decreases until 13 years of rule and rapidly increases after this change point. This change point indicates a change in the Roman Empire’s tendency and dynamics after a given time of reign. Further, we also consider the Byzantine Empire’s time-to-violent-death that, combined with the initial data, comprehends 175 emperors. In this case, we observe a more complex behaviour with additional change points and, therefore, we generalize the power-law model to include multiple change points.

Our model is mainly general and can be used in any system in which the reliability function follows a power law, and we have censored data. Indeed, our work opens the possibility of studying empires and dynasties by using statistical analysis tools, contributing to the understanding of social and historical sciences more quantitatively.

This paper is organized as follows. Initially, we consider the data by Saleh [5] and show that the power-law model is more suitable to model the distribution of time-to-violent-death than previous models. We identify a change point that occurs at 13 years of reign. Next, we consider a dataset of 175 Roman emperors, including also the Byzantine Empire. In this case, we introduce a generalized power-law distribution with k − 1 change points and show that this distribution accurately describes the data observed. Finally, we analyse some possible causes of the short lives of Roman emperors and verify that the birthright is a statistically significant variable, whereas other features, such as the birth province and dynasty, do not influence the time-to-violent-death. These findings suggest that emperors who inherited their reign tend to have a more peaceful reign and increased probability of natural death. At the end of the paper, we present the mathematical formulation and simulation study related to our statistical model.

Survival analysis

The power-law distribution can be represented by its survival function (also refereed by complementary cumulative density function) given bywhere α > 1 and tmin > 0 is the minimum value at which power-law behaviour holds. The inferential methods to estimate α have been discussed by many authors [14-16]. The most common approach to this end is the MLE, which in the continuous case has a simple closed-form expression as well as desired properties, such as consistency and asymptotically efficiency.

To construct our model, let us consider figure 1. The top panel presents the Kaplan–Meier non-parametric estimator for the time-to-violent-death of the 69 rulers of the Roman Empire [5]. The Kaplan–Meier estimator shows the probability of an event at a certain time interval. It can be noticed that the survival time follows a power-law distribution, and after a change point, the behaviour follows a power-law distribution with a different scaling exponent. In figure 1, we consider a power-law distribution with α1 parameter and after a critical value t, the behaviour change to a power-law distribution with coefficient α2.

Figure 1

Illustration of the survival function modelled by a power-law distribution with change point t. The Kaplan–Meier estimator is shown inside the figure.

In a general form, the power-law model behaviour in figure 1 can be represented by the following expression:where is an identity function that returns 1 if the value is inside the interval and 0 otherwise, is the normalized constant that make the model a complementary cumulative function and t > tmin. From equation (2.1) the probability distribution of the power-law distribution with change point is given by

An important aspect of the proposed model is that it can be seen as a generalization of the power-law distribution in the thermodynamic limit of t, i.e.which leads to the standard power-law distribution.

In many cases, it is useful to present other mathematical functions related to the proposed model, such as the mean, variance and kurtosis, to name a few. These functions can be obtained from the rth moment that is derived in appendix A.1. Another aspect that needs to be discussed is the inferential procedure used to obtain the parameter estimates. Although the MLE has already been discussed for the power-law distribution (see Clauset et al. [14] and references therein), this method needs to be generalized to our model. The closed-form estimators for the parameters and the asymptotic confidence intervals obtained from the maximum-likelihood estimators are discussed in detail in appendix A.2.

The estimators presented in appendix A for complete data can be used satisfactorily in most scenarios. However, these estimators may return biased estimates in the presence of events with incomplete information. For instance, Titus rules for less than 1 year and died, probably owing to natural causes, however, he never experienced the event of interest, i.e. ‘a violent death’, in this case, his information is considered as censored, he may or may not experience such type of death if he had ruled for a long period such as Augustus. To take into account such information, we proposed a model that accommodates the presence of two distinct groups: one that suffered a violent death and the other that would not suffer a violent death even if we wait for a long time. These models are known as long-term survival distributions and the parameter estimation that takes into account this characteristic is discussed in detail in appendix A.2.2. This procedure is used to achieve the estimates for the datasets discussed above.

Results

We divide our research into two parts. Initially, we analyse the data by Saleh [5], which considers the western Roman Empire, from Augustus (d. 14 CE) to Theodosius (d. 395 CE). In this analysis, we compare our model with the previous model discussed in [5]. Our study includes the eastern Roman Empire, which ended with Constantine XI Palaiologos (d. 1453 CE). While the first dataset considers only 69 emperors, the second one is complete, covering 175 Roman emperors.

Western Roman Empire

Initially, we consider the data organized by Saleh [5], which was obtained from De Imperatoribus Romanis, an online encyclopedia of Roman emperors available as supplementary material in the cited paper. Here, we consider the same dataset in order to compare the two models. Saleh [5] has considered a mixture of Weibull distributions to describe the lifetime, given by

An important characteristic that should be taken into account when analysing data from the Roman Empire is the existence of emperors who did not have a violent death, i.e. a portion of the emperors was not susceptible to the event of interest. From the model (equation (3.1)), we have that S(t) → 0 as the time increases. Hence, the model assumes that all the emperors will experience a violent death at some time. This assumption is not plausible, as many emperors died owing to natural causes and did not have a violent death. Such characteristics can be incorporated into our power-law model by considering a long-term survival model. In this case, S(t) → π as the time increase, where π ≥ 0 is the proportion of emperors that did not experience a violent death. The details related to the modifications in the power-law distribution with a change point to include long-term survival can be seen in appendix A.4.

Fitting the power-law distribution to the data, the MLEs are obtained and return the estimates , and π = 0.22 when fixing tmin = 0.5 and t = 13. Although these results returned good shapes for the survival function we observe that there is space for improvement. Considering a simple grid search around the obtained values we achieve the following estimates , and π = 0.245 which provided better curves when compared with the empirical survival function. Hence, the adjusted model is given by

Owing to discontinuity in the PDF, the standard maximum likelihood of t cannot be achieved. This problem can be easily overcome by considering a grid search for the parameter where the Kolmogorov–Smirnov test is used to estimate t. Figure 2 presents the fit of our proposed model compared with the mixture Weibull and the non-parametric Kaplan–Meier estimator. Additionally, we constructed a goodness-of-fit test for our proposed model (described in detail in appendix A.6) by using Geerdens et al. [17] method that uses Cramér–von Mises distance and a bootstrap technique to obtain the p-value; for the proposed dataset we obtained a pvalue = 0.99, which indicates that at 5% confidence level our model can be assumed to describe the data.

Figure 2

Power-law survivor (reliability) function of Roman emperors, the adjusted mixture Weibull and the nonparametric Kaplan–Meier estimates. We consider tmin = 0.5, α1 = 1.382, t = 13 and α2 = 8.5 in the power-law model shown in equation (2.1).

Our results indicate that the change point occurs at 13 years of reign. An interesting fact is that, although Alexander the Great (d. 323 BC) has not participated in the unified Roman Empire and was not included in our analysis, his reign had 13 years. From the adjusted results, we observe that Gaius Caligula (d. 41) had a violent death before the change point, while Nero (d. 68 AC) experienced this type of death after the change point. Although we have observed about 37% of censored information, our model estimates that nearly 25% (π = 0.248) of the emperors would not die owing to a violent death. For example, Marcus Aurelius (d. 160 AC), known as the last of the five good emperors, lived in an age of relative peace and died owing to natural causes. As mentioned earlier, Augustus (d. 14) also reigned for 40 years and did not have a violent death. While our model can capture this effect in the survival function, the mixture Weibull distribution cannot return very far estimates when compared with the non-parametric survival, which is undesirable.

The breakpoint was observed at nearly 13 years and may have many different causes about which some speculation can be done. Saleh [5] argued that this period may be related to the mismatch between capabilities and demands under changing political and geopolitical circumstances. Christian & Elbourneb [18] discussed that low periods of rainfall led to increases in the probability that Roman troops starve, making soldiers mutiny and increasing the probability of the emperor being assassinated. Cuaresma et al. analysed the relationship between oil endowment and the duration of dictatorial regimes by comprising information for 106 dictators from 1875 until 2004, they observed a pattern where on average, they were in power for about 12.3 years (closer to our changing point) where 11.3% failed within the first year of their dictatorship; although the time period and the political aspects have changed significantly a similar pattern is observed in dictatorial regimes.

Figure 3 presents the failure rate of the adjusted model using hpop(t) = fpop(t)Spop(t)−1 (given in equations (A 9) and (A 11)) which shows a decreasing failure rate in the first years up to 13 years. During the change point, there is a rapid increase in the hazard of death (violent) and then a decrease again over time. These results differ from the ones obtained from the failure rate of the mixture Weibull distribution. In the latter, the increase in the hazard starts at nearly 11 years, reaches a peak at 15 years, and then decreases. Hence, our proposed approach returns a better understanding of the stochastic process that describes the time-to-violent-death of the emperors of the unified Roman Empire.

Figure 3

Failure rate from power-law with phase transition distribution of Roman emperors. Cox-Snell residuals analysis for the mixture Weibull and power-law distribution.

western and eastern Roman Empire

Although the analysis by Saleh [5] considered the dynasties that started from Augustus (d. 14 CE) and ended with Theodosius (d. 395 CE) when the western Roman Empire began to disintegrate, here, we extend his study considering the period of the Byzantine Empire (eastern Roman Empire) that ended with Constantine XI Palaiologos (d. 1453 CE). The additional data included 106 emperors, where the observed censorship was 56.57%, which is higher than the observed in the first analysed period, and the long-term survival was estimated in π = 0.404, i.e. 40.40% of the Roman emperors probably would not face a violent death. This finding implies that during the Byzantine Empire, the emperors were subjected to a less violent death. This change in the hazard behaviour may be explained owing to different causes; for instance, the eastern half of the Roman Empire was shown to be less susceptible to external attacks owing to its geography. Additionally, the Empire did not have many common boundaries with Europe. Although many other complex causes are also responsible for providing such stability, those facts combined with the internal political cohesion and more robust administrative centre allowed the Byzantine Empire to survive for many centuries.

Here, the survival function related to the time-to-violent-death for all emperors has more complex behaviour, therefore, we have generalized the power-law distribution with k − 1 change points that can be represented by the following expression:where is an identity function that returns 1 if the value is inside the interval and 0 otherwise, C0 = 1, t*( = ∞, α > 0 and . From equation (3.2) the probability distribution of the piecewise power-law model distribution is given by

The details related to the generalized model with long-term survival and the inference procedure to estimate the parameters can be seen in appendix A. This generalization is very flexible, and the fit improves following the number of change points. The model is essentially a mixture of truncated Pareto distributions [19], although some forms of mixture have been presented earlier [20,21] they differ from our proposed model, which is not constructed based on normalized truncated Pareto distributions, in this case, the proposed approach can be seen as a piecewise power-law distribution. Additionally, there is significant literature on change point models [22-24], however, most of the results are focused on the presence of one change point in the hazard function, while other models are based on the assumption that the hazard function follows a piecewise exponential model [25,26], where the hazard rate inside each change-point are assumed to be constant, which is not the case for the proposed data. Here our aim is to show that the proposed data can be described by power-law distributions which is an important pattern that occurs in many complex systems. Figure 4 presents the survival of the adjusted model using k = 4, 5 as well as the presence of long-term survival.

Figure 4

Power-law survivor (reliability) function of Roman emperors for all 175 Roman emperors and the non-parametric Kaplan–Meier estimates.

As shown in figure 4, as the number of change points increases, we obtain a better fit for the proposed data. In fact, with k = 5, the adjusted model almost overlaps the non-parametric function. In order to avoid overfitting, we consider the Akaike information criteria (AIC) to discriminate between the two models. In this case, the best fit among the chosen models returns the minimum value of AIC. While the model with k = 4 returned an AIC = 594.10, the model with k = 5 yields AIC = 586.25. Thus, this is the best model for describing the time-to-violent-death for all emperors. The goodness-of-fit test was also applied for the model which for the extended dataset we obtained a pvalue = 0.337, which indicates that at 5% confidence level our model can be assumed to describe the data.

The time-to-violent-death of western emperors

We also analyse the influence of some emperors’ attributes on the time-to-violent-death. These attributes are compiled from Wikipedia and previous works, including [3,27], covering western emperors from 26 BC to 395 AD. In the case of the method of accession to power (e.g. birthright, seized power, appointment by the senate, among others), we can see in figure 5 that emperors who came from birthright tend to rule longer than from other methods. To verify the significance of this difference, we consider a two-sided test for the null hypothesis that the average time-to-violent-death is independent of the method of accession to power. If the p-value is smaller than the threshold, e.g. 1% or 5%, then we reject the null hypothesis of equal averages. From this hypothesis test, we conclude that the difference between the time-to-violent-death of emperors that used different methods of accession to power is significant because p < 0.001. This result suggests that emperors who inherited the reign tend to have a more peaceful administration and an increased probability of having a natural death. On the other hand, in terms of the birth province, there is no difference between the time-to-violent-death of emperors from Italy or other parts of the world, like Syria or Spain (p = 0.54). This finding indicates that Rome was very cosmopolitan, as the Empire was extensive, including many different cultures, traditions and faiths. Thus, it is expected that the place of birth plays a small part in the time-to-violent-death.

Figure 5

Box plot of emperors’ attributes according to the time-to-death.

The cause of death also does not influence the time-to-violent-death. Indeed, emperors, who died owing to assassination, did not have shorter ruling periods (p = 0.8), although almost 40% of the emperors were assassinated. Related to this, emperors killed by other emperors also did not have a shorter reign period (p = 0.4), as shown in the figure.

In terms of era, we verify the presence of a moderate difference (p = 0.08) in terms of the reign period between the Principate. That is, the first period of the Roman Empire from the beginning of Augustus’s reign in 27 BC, and Dominate, which started after the Crisis of the Third Century, in AD 284. The reigning period is more extended during the Dominate era, with an average time of ruling of 11 years, contrasting with 7 years, verified in the Principate period covered in our dataset.

Regarding the dynasty, we also do not observe a significant difference between them in terms of the reign period (p = 0.2). In figure 5, we show the comparison between the Severan dynasty, which started with Septimius Severus, who rose to power as the victor of the Civil War of 193–197, and other dynasties. Although the Severan dynasty was disturbed by highly unstable family relationships and constant political turmoil [27], it did not significantly affect the time-to-violent-death.

Another important tool to analyse such data is using the Cox proportional hazards model [28], in this case, we can verify if the attributes affect the risk of death in a multiplicative way or provide higher risks during a specific time. The proportional assumption is checked using a test based on Schoenfeld residuals [29] where we present the p-values, assuming the threshold of 5% there is an indication of a violation of the proportional hazard assumption when the values are smaller than 0.05. The analysed attributes and its respective p-values in parenthesis are: birthright (0.2371), birth province (0.5000), cause (0.8785), dynasty (0.7429), era (0.2452), killer (0.0386), global (0.2554). Hence, the p-values indicate that only the attribute killer was shown to influence the change in the risk of death, while other covariates were shown to have proportional hazard rate differing only in a multiplicative term in the hazard function. The emperors who died by other emperors showed a higher risk of death in the beginning of their reign, on the other hand, after 16 years the curves of risks crossed each other, exhibiting a higher risk of dying by others instead of emperors. Overall, most of the attributes considered in our analysis do not affect the time-to-violent-death, which suggests that other features, like those related to external factors like wars and political conflicts, should be considered to infer the causes of the short reigning period of Roman emperors. This understanding is an exciting topic for further research.

Discussion

In this paper, we introduce a new survival model with k change points and power-law distribution. The model is a modification of the traditional power-law model, in which we have included the adjustments to address the long-term survival and the new procedures to obtain the parameter estimates. The model has been developed to predict the reign period of Roman emperors, motivated by the analysis provided by Saleh [5].

In examining the time-to-violent-death of the emperors of the Roman Empire, we have observed that the risk is high when the emperor assumes the throne. This finding may be related to the struggles in handling the demands that the position requires and the lack of political skills of the new emperor. Our model suggests that the risk systematically decreases until 13 years of rule and then rapidly increases after this change point. The are many possible reasons for the change in the behaviour of the risk, for instance, their old enemies had regrouped, or new ones emerged [5]. After the change point, the risk decreases again. Furthermore, the adjusted survival functions lead to 25% of emperors that would not die owing to violent death. Therefore, being a Roman emperor was a hazardous occupation that led to three out of four suffering from a cruel death during the western Empire. Additionally, taking into consideration the eastern Roman Empire, we observed that out of three emperors suffer from a cruel death.

The comparison with the model by Saleh [5] shows that the power-law model is more accurate to describe the variable of interest, i.e. the time-to-violent-death of the emperors of the Roman Empire. This fact has important implications since the power-law distribution is ubiquitous in natural and artificial systems and is often related to critical events. Power laws can be observed in earthquakes [7-9], in solar radiophysics [10], stock price fluctuations [11] and tree-limb branching [12]. In the dataset considered here, the lifetime of Roman emperors decay as a power law, which implies that short reigns are common, while high occurrences are rare. This behaviour is similar to that observed in many systems, including earthquake, where huge earthquakes are rare. This result suggests that the time-to-violent-death does not have a typical scale, including the Roman Empire as an example of systems presenting scale invariances, such as the distribution of income, size of cities according to population, size of corporations and word frequencies [14].

Our power-law survival model with k change points is general and can be applied to model the hazard rate function in engineering, medicine, economy, and physics, especially in situations where we have a good indication that the tail of the distribution follows a power law. However, the earlier occurrences may have a more slow decay. In addition to the proposed results drawn here, there are additional aspects of the new distribution that should be considered. For instance, the theoretical results assuming discrete data are currently under investigation. Under this scenario, the estimation procedures do not have a closed-form expression and need to be further improved. Thus, our results open several possibilities in applications and theoretical development, including models with cut-off and improvements by considering Bayesian methods.

Click here for additional data file.