Explanation of machine learning models is an important problem in many applications. For instance, medicine machine learning applications meet a requirement of understanding results provided by the models. A typical example is when a doctor has to have an explanation of a stated diagnosis in order to have an opportunity to choose a preferable treatment (Holzinger et al., 2019). This implies that decisions provided by the machine learning models should be explainable to help machine learning users to understand the obtained results, for example, doctors need to understand obtained diagnoses. One of the obstacles to obtain explainable decisions is the black-box nature of many models, especially, of deep learning models, i.e. inputs and outcomes of these models may be known, but there is no information what features impact on corresponding decisions provided by the models. Many explanation methods have been developed in order to overcome this obstacle and to explain the model outcomes. The explanation methods can be divided into two groups: local and global methods. Methods from the first group derive explanation locally around a test example, whereas methods from the second group try to explain the black-box model on the whole dataset or a part of the datasets. The global explanation methods are of great interest, but many application areas, especially, medicine, require to understand decisions concerning a specific patient, i.e. it is important to know what features of an example are responsible for a black-box model prediction. Therefore, this paper focuses on local explanations.

It is important to note that users of a black-box model are often not interested why a certain prediction was obtained and what features of an example led to a decision. They usually aim to understand what could be changed to get a preferable result by using the black-box model. Such explanations can be referred to the so-called counterfactual explanations or counterfactuals (Wachter et al., 2017), which may be more desirable, intuitive and useful than “direct” explanations (methods based on attributing a prediction to input features). According to Molnar (2019), a counterfactual explanation of a prediction can be defined as the smallest change to the feature values of an input original example that changes the prediction to a predefined outcome. There is a classic example of the loan application rejection (Molnar, 2019; Wachter et al., 2017), which explicitly explains counterfactuals. According to this example, a bank rejects an application of a user for a credit. A counterfactual explanation could be that “the credit would have been approved if the user would earn $500 more per month and have the credit score 30 points higher” (Molnar, 2019; Wachter et al., 2017).

So far, methods of counterfactual explanations have been applied to classification and regression problems where a black-box model produces a point-valued outcome for every input example. However, there are many models where the outcome is a function. Some of these models solve problems of survival analysis (Hosmer et al., 2008), where the outcome is a survival function (SF) or a cumulative hazard function (CHF). In contrast to the standard classification and regression machine learning models, the survival models deal with datasets containing a lot of censored data. This fact complicates the models.

This paper presents a method for finding counterfactual explanations for predictions of survival machine learning black-box models, which is based on analysis of SFs. The method allows us to find a counterfactual whose SF is connected with the SF of the original example by means of some conditions. One of the difficulties of solving the counterfactual explanation problem is that the classes of examples are implicitly defined through outcomes of a machine learning survival model in the form of SFs or CHFs. Therefore, a condition establishing the difference between mean times to event of the original example and the counterfactual is proposed. For example, the mean time to event of the counterfactual should be larger than the mean time to event of the original example by value r. The meaning of counterfactuals in survival analysis under the above condition can be represented by the statement: “Your treatment was not successful because of a small dose of the medicine (one tablet). If your dose had been three tablets, the mean time of recession would have been increased to a required value”. Here the difference between the required value of the mean time to recession and the recent mean time to recession of the patient is r. It depends on the black-box model used in a certain study. In particular, when the Cox model is considered as a black-box model, the exact solution can be obtained because the optimization problem for computing counterfactuals is reduced to a standard convex programming problem. In a general case of the black-box model, the optimization problem for computing counterfactuals is non-convex. Therefore, one of the ways for solving the optimization problem is to use some heuristic global optimization method. An optimization method called Particle Swarm Optimization (PSO), proposed by Kennedy and Eberhart (1995), can be applied to solving the counterfactual explanation problem. The method is a population-based stochastic optimization technique based on swarms and motivated by the intelligent collective behavior of some animals (Wang et al., 2018).

In summary, the following contributions are made in this paper:

The paper is organized as follows. Related work is briefly discussed in Section 2. Basic concepts of survival analysis and the Cox model are given in Section 3. Section 4 contains the standard counterfactual explanation problem statement and its extension to the case of survival analysis. In Section 5, it is shown that the counterfactual explanation problem for the black-box Cox model is a convex programming problem and therefore can simply be solved. Its application to the counterfactual explanation problem is considered in Section 7. Numerical experiments with synthetic data and real data are given in Section 8. Discussion of some peculiarities of the proposed method and concluding remarks can be found in Section 9.

Local explanation methods and counterfactual explanations. Due to importance of the machine learning model explanation in many applications, many methods have been proposed to explain black-box models locally (Arya et al., 2019; Guidotti et al., 2019b; Molnar, 2019; Murdoch et al., 2019). A critical review and analysis of many explanation methods can be found in survey papers (Adadi and Berrada, 2018; Arrieta et al., 2019; Carvalho et al., 2019; Das and Rad, 2020; Guidotti et al., 2019b; Rudin, 2019; Xie et al., 2020). One of the first local explanation methods is the LIME method (Ribeiro et al., 2016; Garreau and von Luxburg, 2020), which uses simple and easily understandable linear models to approximate the predictions of black-box models locally. LIME, as well as many local explanation methods, are based on perturbation techniques (Fong and Vedaldi, 2019, 2017; Petsiuk et al., 2018; Vu et al., 2019). Another explanation method, which is based on the linear approximation, is the SHAP (Lundberg and Lee, 2017; Strumbelj and Kononenko, 2010), which takes a game-theoretic approach for optimizing a regression loss function based on Shapley values.

In order to get intuitive and human-friendly explanations, several counterfactual explanation methods (Wachter et al., 2017) were developed by several authors (Buhrmester et al., 2019; Dandl et al., 2020; Fernandez et al., 2020; Goyal et al., 2019; Guidotti et al., 2019a; Hendricks et al., 2018b; Looveren and Klaise, 2019; Lucic et al., 2019; Poyiadzi et al., 2020; Russel, 2019; Vermeire and Martens, 2020; van der Waa et al., 2018; White and Garcez, 2019). The counterfactual explanation tells us what to do to achieve a desired outcome.

Some counterfactual explanation methods use combinations with other methods like LIME and SHAP. For example, (Ramon et al., 2019) proposed the so-called LIME-C and SHAP-C methods as counterfactual extensions of LIME and SHAP. White and Garcez (2019) proposed the CLEAR methods which can also be viewed as a combination of LIME and counterfactual explanations.

Many counterfactual explanation methods apply perturbation techniques to examine feature perturbations that lead to a different outcome of a black-box model. In fact, this is an approach to generate counterfactuals to alter an input of the black-box model and to observe how the output changes. One of the methods using perturbations is the Growing Spheres method (Laugel et al., 2018), which can be referred to counterfactual explanations to some extent. The method determines the minimal changes needed to alter a prediction. Perturbations are usually performed towards interpretable counterfactuals in a lot of methods (Dhurandhar et al., 2018, 2019; Looveren and Klaise, 2019).

Another direction for applying counterfactuals concerns with counterfactual visual explanations which can be generated to explain the decisions of a deep learning system by identifying what and how regions of an input image would need to change in order for the system to produce a specified output (Goyal et al., 2019). Hendricks et al. (2018a) proposed counterfactual explanations in natural language by generating counterfactual textual evidence. Counterfactual “feature-highlighting explanations” were proposed by Barocas et al. (2020) by highlighting a set of features deemed most relevant and withholding others. A counterfactual impact analysis of medical images was considered by Lenis et al. (2020), and by Bhatt et al. (2019).

Many other approaches have been proposed in the last few years (Verma et al., 2020), but there are no counterfactual explanations of predictions provided by the survival machine learning systems. Therefore, our aim is to propose a new method for counterfactual survival explanations.

Machine learning models in survival analysis. Survival analysis is an important direction for taking into account censored data. It covers many real application problems, especially in medicine, reliability analysis, risk analysis. Therefore, the machine learning approach for solving survival analysis problems allows improving the survival data processing. Many machine learning survival models have been developed (Lee et al., 2018; Wang et al., 2019; Zhao and Feng, 2020) due to importance of survival models in many applications, including reliability of complex systems, medicine, risk analysis, etc. A comprehensive review of machine learning models dealing with survival analysis problems was provided by Wang et al. (2019). One of the most powerful and popular methods for dealing with survival data is the semi-parametric Cox proportional hazards model (Cox, 1972). Many modifications have been developed to relax some strong assumptions underlying the Cox model. In order to take into account the high dimensionality of survival data and to solve the feature selection problem with these data, Tibshirani (1997) presented a modification based on the Lasso method. Similar Lasso modifications, for example, the adaptive Lasso, were also proposed by several authors (Kim et al., 2012; Witten and Tibshirani, 2010; Zhang and Lu, 2007). A further extension of the Cox model is a set of SVM modifications (Khan and Zubek, 2008; Widodo and Yang, 2011). Various architectures of neural networks, starting from a simple network (Faraggi and Simon, 1995), proposed to relax the linear relationship assumption in the Cox model, have been developed (Haarburger et al., 2018; Katzman et al., 2018; Ranganath et al., 2016; Zhu et al., 2016) to solve prediction problems in the framework of survival analysis. Despite many powerful machine learning approaches for solving the survival problems, the most efficient and popular tool for survival analysis under condition of small survival data is the extension of the standard random forest (Breiman, 2001) called the random survival forest (RSF) (Ibrahim et al., 2008; Mogensen et al., 2012; Wang and Zhou, 2017; Wright et al., 2017).

Most of the above models dealing with survival data can be regarded as black-box models and should be explained. However, only the Cox model has a simple explanation due to its linear relationship between covariates. Therefore, it can be used to approximate more powerful models, including survival deep neural networks and random survival forests, for explaining predictions of these models.

Explanation methods in survival analysis. There are several methods explaining survival machine learning model predictions. Kovalev et al. (2020) proposed an explanation method called SurvLIME, which deals with censored data and can be regarded as an extension of the Local Interpretable Model-agnostic Explanations (LIME) method (Ribeiro et al., 2016) on the case of survival data. The basic idea behind SurvLIME is to apply the Cox model to approximate the black-box survival model at a local area around a test example. SurvLIME used the quadratic norm to take into account the distance between CHFs. Following ideas behind SurvLIME, Utkin et al. (2020) proposed a modification of SurvLIME called SurvLIME-Inf. In contrast to SurvLIME, SurvLIME-Inf uses L ∞ -norm for defining distances between CHFs. SurvLIME-Inf significantly simplifies the model and provides better results when a training set is small. Another explanation model proposed by Kovalev and Utkin (2020) is called SurvLIME-KS. This model uses the well-known Kolmogorov-Smirnov bounds to ensure robustness of the explanation model in cases with a small amount of training data or outliers of survival data.

We consider a definition of the counterfactual explanation proposed by Wachter et al. (2017) and rewrite it in terms of survival models.

The function l ( f ( z ) , f ( x ) ) encourages the prediction of z to be different in accordance with a certain rule to the prediction of the original point x. The penalty term μ ( z , x ) minimizes the distance between z and x with the aim to find the nearest counterfactuals to x. It can be defined as (8) μ ( z , x ) = ‖ z − x ‖ 2 .

It is important to note that the above optimization problem can be extended by including additional terms. In particular, many algorithms of the counterfactual explanation use a term which makes counterfactuals close to the observed data. It can be done, for example, by minimizing the distance between the counterfactual z and the k nearest observed data points (Dandl et al., 2020) or by minimizing the distance between the counterfactual z and the class prototypes (Looveren and Klaise, 2019).

Let us consider an analogy of survival models with the standard classification models where all points are divided into classes. We also have to divide all patients into classes by means of an implicit relationship between the black-box survival model predictions. It is important to note that predictions are the CHFs or the SFs. Therefore, the introduced loss function l ( f ( z ) , f ( x ) ) should take into account the difference between the CHFs or the SFs to some extent, which characterize different “classes” or groups of patients. It is necessary to establish the relationship between CHFs or between SFs, which would separate groups of patients of interest. One of the simplest ways is to separate groups of patients in accordance with their feature vectors. This can be done if the groups of patients are known, for example, the treatment and control groups. In many cases, it is difficult to divide patients into groups by relying on their features because this division does not take into account outcomes, for example, SFs of patients.

Another way for separating patients is to consider the difference between the corresponding mean times to events for counterfactual z and input x. Therefore, several conditions of counterfactuals taking into account mean values can be proposed. The mean values can be defined as follows: (9) m ( z ) = E ( z ) = ∫ 0 ∞ S ( t | z ) d t , m ( x ) = E ( x ) = ∫ 0 ∞ S ( t | x ) d t .

Then the optimization problem (7) can be rewritten as follows: (10) min z ∈ R d { l ( m ( z ) , m ( x ) ) + C μ ( z , x ) } .

We suppose that a condition for a boundary of “classes” of patients can be defined by a predefined smallest distance between mean values, which is equal to r. In other words, a counterfactual z is defined by the following condition: (11) m ( z ) − m ( x ) ⩾ r .

The condition for “classes” of patients can be also written as (12) m ( x ) − m ( z ) ⩾ r .

Let us unite the above conditions by means of the function (13) ψ ( z ) = r − θ ( m ( z ) − m ( x ) ) ⩽ 0 , where parameter θ ∈ { − 1 , 1 }. In particular, condition (11) corresponds to case θ = 1, condition (12) corresponds to case θ = − 1.

It should be noted that several conditions of counterfactuals taking into account mean values can be proposed here. We take the difference between the mean time to events of explainable point x and the point z. For example, if it is known that a group of patients has a certain disease, we can define a prototype x p of the group of patients and try to find the difference m ( z ) − m ( x p ).

Let us consider the so-called hinge-loss function (14) l ( f ( z ) , f ( x ) ) = max ( 0 , r − ( f ( z ) − f ( x ) ) ) .

Its minimization encourages to increase the difference f ( z ) − f ( x ) up to r. Indeed, the condition r ⩾ 0 is valid, therefore, the increase of f ( z ) leads to the decrease of l ( f ( z ) , f ( x ) ). However, when f ( z ) − f ( x ) ⩾ r, then l ( f ( z ) , f ( x ) ) = 0. This implies that minimization of l ( f ( z ) , f ( x ) ) does not encourage to increase the corresponding difference. But it does not mean that f ( z ) − f ( x ) cannot be larger than r. In fact, the loss function minimization encourages to move point f ( z ) to the class boundary, but it does not impact on its moving in the area of another class.

Taking into account (8) and (13), the entire loss function can be rewritten and the following optimization problem is formulated: (15) min z ∈ R m L ( z ) = min z ∈ R m { max { 0 , ψ ( z ) } + C ‖ z − x ‖ 2 } .

In summary, the optimization problem (15) has to be solved in order to find counterfactuals z.

It should be noted that counterfactuals can also be found by solving the following constrained optimization problem: (16) min z ∈ R m L ( z ) = min z ∈ R m ‖ z − x ‖ 2 , subject to (17) ψ ( z ) ⩽ 0 .

It is equivalent to problem (15) or to problem (7), which are simply derived from (16)–(17).

Generally, the function ψ ( z ) is not convex and cannot be written in an explicit form. This fact complicates the problem and restricts possible methods for its solution. Therefore, we propose to use the well-known heuristic method called the PSO. Let us return to the problem (15) by writing it in the similar form: (18) min z ∈ R m L ( z ) = min z ∈ R m { max { 0 , C ψ ( z ) } + ‖ z − x ‖ 2 } .

This problem can be explained as follows. If point z is from the feasible set defined by condition ψ ( z ) ⩽ 0, then the choice of z minimizes the distance between z and x. If point z does not belong to the feasible set ( ψ ( z ) > 00$]]>), then the penalty C ψ ( z ) is assigned. It is assumed that the value of C ψ ( z ) is much larger than ‖ z − x ‖ 2 . Therefore, it is recommended to take large values of C, for example, C = 10 6 .

We prove below that problem (18) can be reduced to a standard convex optimization problem with linear constraints and can be exactly solved if the black-box model is the Cox model. In other words, the counterfactual example z can be determined by solving a convex optimization problem.

Let t 0 < t 1 < ⋯ < t q < t q + 1 be the distinct times to event of interest from the set { T 1 , … , T n }, where t 0 = 0, t 1 = min k = 1 , … , n T k and t q = max k = 1 , … , n T k , t q + 1 = t q + t γ ; t γ is a parameter which is close to 0. We assume that there hold S ( τ | x ) = 1 for τ = t 0 , 0 < S ( τ | x ) < 1 for τ ⩾ t 1 , and S ( τ | x ) = 0 for τ > t q + 1 {t_{q+1}}$]]>. Let Ω = [ t 0 , t q + 1 ] and divide it into q + 1 subsets Ω 0 , … , Ω q such that Ω q = [ t q , t q + 1 ], Ω j = [ t j , t j + 1 ), Ω = ⋃ j = 0 , … , q Ω j ; Ω j ∩ Ω k = ∅, ∀ j ≠ k. Introduce the indicator function χ j ( t ) taking the value 1 if t ∈ Ω j , and 0 otherwise. Then the baseline SF S 0 ( τ ) and the SF S ( τ | x ) under condition of using the Cox model can be represented as follows: (19) S 0 ( τ ) = ∑ j = 0 q s 0 , j · χ j ( τ ) , s 0 , 0 = 1 , and (20) S ( τ | x ) = ∑ j = 0 q s 0 , j exp ( z T b ) · χ j ( τ ) .

Hence, the mean value is (21) m ( x ) = ∫ 0 ∞ S ( t | x ) d t = ∫ 0 ∞ [ ∑ j = 0 q s 0 , j exp ( x T b ) χ j ( t ) ] d t = ∑ j = 0 q s 0 , j exp ( x T b ) [ ∫ 0 ∞ χ j ( t ) d t ] = ∑ j = 0 q μ j s 0 , j exp ( x T b ) , where μ j = t j + 1 − t j > 00$]]>.

Denote u = z T b and consider the function (22) π ( u ) = ∑ j = 0 q μ j s 0 , j exp ( u ) .

Compute the following limits: (23) lim u → − ∞ π ( u ) = ∑ j = 0 q μ j lim u → − ∞ s 0 , j exp ( u ) = ∑ j = 0 q μ j = t q + 1 − t 0 = t q + 1 , (24) lim u → ∞ π ( u ) = ∑ j = 0 q μ j lim u → ∞ s 0 , j exp ( u ) = μ 0 = t 1 − t 0 = t 1 .

The derivative of π ( u ) is (25) d π ( u ) d u = ∑ j = 0 q μ j d d u [ s 0 , j exp ( u ) ] = ∑ j = 0 q [ μ j ln ( s 0 , j ) ] ( s 0 , j exp ( u ) exp ( u ) ) .

Note that s 0 , j exp ( u ) exp ( u ) ⩾ 0 for all j and u; μ j ln ( s 0 , j ) ⩽ 0 for all j. Hence, there holds (26) d π ( u ) d u ⩽ 0 , ∀ u .

The above means that the function π ( u ) is non-increasing with u. Moreover, it is positive because its limits are positive too. Let us consider the function (27) ζ ( u ) = r − θ ( π ( u ) − m ( x ) ) .

It is obvious that there holds m ( x ) ∈ [ t 1 , t q ] for arbitrary x.

Let θ = 1. Then

It follows from the above that a + ⩽ 0 < b + and ζ ( u ) = 0 at a single point u + . By using numerical methods, we can find point u + . Since the set of solutions is defined by the inequality ψ ( z ) ⩽ 0, then it is equivalent to ζ ( u ) ⩽ 0 and u ⩽ u + or z T b − u + ⩽ 0.

Let θ = − 1. Then

It follows from the above that b − ⩽ 0 < a − and ζ ( u ) = 0 at a single point u − which can be numerically calculated. Since the set of solutions is defined by the inequality ψ ( z ) ⩽ 0, then it is equivalent to ζ ( u ) ⩽ 0 and u ⩾ u − or − z T b + u − ⩽ 0.

The above conditions for r and the sets of solutions can be united (28) r ∈ ( 0 , 1 2 [ ( 1 + θ ) ( t q + 1 − m ( x ) ) + ( 1 − θ ) ( m ( x ) − t 1 ) ] ] , (29) θ z T b − 1 2 [ ( 1 + θ ) u + − ( 1 − θ ) u − ] ⩽ 0 .

Constraints (17) to the problem (16)–(17) become (29) which are linear with z. It can be seen from the objective function (16) and constraints (29) that this optimization problem is convex with linear constraints, therefore, it can be solved by means of standard programming methods.

The problem (16) with constraints (29) can also be written in the form of the unconstrained problem (18), as follows: (30) min z ∈ R m L ( z ) = min z ∈ R m { C ( θ z T b − A ) + ‖ z − x ‖ 2 } , where (31) A = 1 2 [ ( 1 + θ ) u + − ( 1 − θ ) u − ] .

The PSO algorithm proposed by Kennedy and Eberhart (1995) can be viewed as a stochastic optimization technique based on a swarm. There are several survey papers devoted to the PSO algorithms, for example, Wang et al. (2018, 2015). We briefly introduce this algorithm below.

The PSO performs searching via a swarm of particles that updates from iteration to iteration. In order to reach the optimal or suboptimal solution to the optimization problem, each particle moves in the direction to its previously best position (denoted as “pbest”) and the global best position (denoted as “gbest”) in the swarm. Suppose that the function f : R n → R has to be minimized. The PSO is implemented in the form of the following algorithm:

The problem has five parameters: the number of particles N; the number of iterations N iter ; the inertia weight w; the coefficient for the cognitive term c 1 (the cognitive term helps the particles for exploring the search space); the coefficient for the social term c 2 (the social term helps the particles for exploiting the search space).

It is clear that parameters N and N iter should be as large as possible. Their upper bounds depend only on the available computation time that we can spend on iterations. We take N = 2000, N iter = 1000.

Other parameters are selected by using some heuristics (Bai, 2010; Clerc and Kennedy, 2002), namely, (36) w = η , c 1 = η ϕ 1 , c 2 = η ϕ 2 , where (37) η = 2 κ | 2 − ϕ − ϕ 2 − 4 ϕ | , ϕ = ϕ 1 + ϕ 2 > 4 , κ ∈ [ 0 , 1 ] . 4,\hspace{2.5pt}\kappa \in [0,1].\]]]>

The following values of the above introduced parameters are often taken: ϕ 1 = ϕ 2 = 2.05, κ = 1. Hence, there hold w = 0.729, c 1 = c 2 = 1.4945.

Particles are generated by means of the uniform distributions U with the following parameters: (38) ( u k 0 ) i ∼ U ( u i min , u i max ) , u i min , u i max ∈ R .

Velocities are similarly generated as: (39) ( v k 0 ) i ∼ U ( − | u i max − u i min | , − | u i max − u i min | ) .

It should be noted that PSO is similar to the Genetic Algorithm (GA) in the sense that they are both population-based search approaches and that they both depend on information sharing among their population members to enhance their search processes using a combination of deterministic and probabilistic rules. However, many authors (Duan et al., 2009; Panda and Padhy, 2008; Sooda and Nair, 2011; Wang et al., 2008) claim that PSO has a better performance in terms of average and standard deviation from multiple runs of algorithms. PSO converges to arrive at the optimal values in fewer generations than GA. Moreover, PSO outperforms GA, when a smaller population size is available, and has higher robustness.

Let us return to the counterfactual explanation problem in the framework of survival analysis. Suppose that there exists a dataset D with triplets ( x j , δ j , T j ), where x j ∈ R d , T j > 00$]]>, δ j ∈ { 0 , 1 }. It is assumed that the explained machine learning model q ( x ) is trained on D. It should be noted that the prediction of the machine learning survival model is the SF or the CHF, which can be used for computing the mean time to event of interest m ( x ). In order to find the counterfactual z, we have to solve the optimization problem (18) with fixed x, r, and C.

Let us calculate bounds of the domain x for every feature on the basis of the training set as (40) x i min = min j { ( x j ) i } , x i max = max j { ( x j ) i } .

According to the PSO algorithm, the initial positions of particles are generated as (41) ( u k 0 ) i ∼ U ( x i min , x i max ) .

So, the optimal solution can be found in the hyper parallelepiped X: (42) X = [ x 1 min , x 1 max ] × ⋯ × [ x d min , x d max ] .

If there exists at least one point x j ∗ in the training set such that ψ ( x j ∗ ) ⩽ 0, then the region X can be adjusted. Let (43) z closest , train = z c t = arg min j L ( x j ) .

Let us introduce a sphere B = B ( x , R c t ) with centre x and radius R c t = ‖ x − z c t ‖ 2 . The sphere can be partially located inside the hyper parallelepiped X or can be larger. Therefore, we restrict the set of solutions by a set M defined as M = X ∩ B. To disable a possible passage beyond the limits of M, we introduce the restriction procedure, denoted as RP, which supports that:

The initial positions of particles are generated as follows: (44) u 1 0 = z c t , u 2 0 , … , u N 0 ∼ U ( B ) , and the above restriction procedure is used for all points except the first one: u k 0 = RP ( u k 0 ), k = 2 , … , N. Positions of particles will be adjusted by using the following expression: (45) u k t = RP ( u k t − 1 + v k t ) .

Initial values of velocities are taken as v k 0 = 0, k = 1 , … , N.

Let us point out properties of the above approach:

Another important question arising with respect to the above approach on the basis of the PSO is how to take into account categorical features. We have to note that the proposed method can potentially deal with categorical features. A direct way for taking into account these features is to consider the optimization problem (18) for different combinations of values of categorical features. Let us represent the feature vector z as ( z ∖ z c a t , z c a t ), where z c a t is the vector consisting of c categorical features. This implies that the problem (18) has to be solved d 1 · d 2 · … · d c times where d i is the number of values of the i-th categorical feature. The counterfactual can be found by minimizing the loss function min z ∖ z c a t L ( z | z c a t ) over all possible vectors z c a t , i.e. we have to find a vector z c a t of the categorical feature values, which provides the smallest value of min z ∖ z c a t L ( z | z c a t ). Here, L ( z | z c a t ) is the loss function under condition of fixed values of z c a t . The above approach is obvious and can be applied to finding the counterfactual in the case of a small number of categorical features.

To perform numerical experiments, we use the following general scheme.

1. The Cox model and the RSF are considered as black-box models that are trained on synthetic or real survival data. Outputs of the trained models in the testing phase are SFs.

2. To study the proposed explanation algorithm by means of synthetic data, we generate random survival times to events by using the Cox model estimates.

In order to analyse the numerical results, the following schemes are proposed for verification. When the Cox model is used as a black-box model, we can get the exact solution. This implies that we can exactly compute the counterfactual z ver . Adding the condition that the solution belongs to the hyper parallelepiped X to the problem with objective function (16) and constraints (29), we use this solution ( z ver ) as a referenced solution in order to compare another solution ( z opt ) obtained by means of the PSO. The Euclidean distance between z ver and z opt can be a measure for the PSO algorithm accuracy in the case of the black-box Cox model.

The next question is how to verify results of the RSF as a black-box model. The problem is that the RSF does not allow us to obtain exact results by means of formal methods, for example, by solving the optimization problem (16). However, the counterfactual can be found with arbitrary accuracy by considering all points or many points in accordance with a grid. Then the minimal distance between the original point x and each generated point is minimized under condition ψ ( z ) ⩽ 0 which is verified for every generated z. This is a computationally expensive task, but it can be applied to testing results. By using the above approach, many random points are generated from the set M defined in the previous section and approximate the optimum z ver . Random points for verification of results obtained by using the RSF are uniformly selected from sphere B by using the restriction procedure RP. The number of the points is set at 10 6 . In fact, this approach can be regarded as a perturbation method with the exhaustive search. The Euclidean distance between z ver and z opt is the accuracy measure when the black-box model is the RSF, but z ver in this case has another meaning than in the case of the black-box Cox model.

The code of the proposed algorithm in Python is available at https://github.com/kovmax/XAI_Survival_Counterfactual.