Counterfactual explanation of machine learning survival models

A method for counterfactual explanation of machine learning survival models is proposed. 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 survival functions. A condition that establishes the difference between survival functions of the original example and the counterfactual is introduced. This condition is based on using a distance between mean times to event. It is shown that the counterfactual explanation problem can be reduced to a standard convex optimization problem with linear constraints when the explained black-box model is the Cox model. For other black-box models, it is proposed to apply the well-known Particle Swarm Optimization algorithm. A lot of numerical experiments with real and synthetic data demonstrate the proposed method.


Introduction
Explanation of machine learning models has become an important problem in many applications. For instance, a lot of 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 [29]. This implies that the machine learning model decisions can be trusted and explainable to help machine learning users or experts to understand the obtained results. 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. A lot of explanation methods have been developed in order to overcome this obstacle [4,25,46,47] 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 its part. The global explanation methods are of a high interest, but many application areas, especially, medicine, require to understand decisions concerning with a 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 [66], which may be more desirable, intuitive and useful than "direct" explanations (methods based on attributing a prediction to input features). According to [46], 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 [46,66], 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" [46,66].
So far, methods of counterfactual explanations have been applied to classification and regression problems where a black-box model produces a point-value outcome for every input example. However, there are a lot of models, where the outcome is a function. A part of these models solves problems of survival analysis [30], where the outcome is a survival function (SF) or a cumulative hazard function (CHF). In contrast to usual classification and regression machine learning models, the survival models deal with datasets containing a lot of censored data. This fact complicates the models. A lot of machine learning survival models have been developed [41,69,77] due to importance the survival models in many applications, including reliability of complex systems, medicine, risk analysis, etc. In spite of some strong assumptions, the semi-parametric Cox proportional hazards model [13] remains one of the most popular one. It is one of the models establishing an explicit relationship between the covariates and the distribution of survival times. The Cox model assumes a linear combination of the example covariates. On the one hand, this is a too strong assumption that is not valid in many cases. It restricts the wide use of the model. On the other hand, this assumption allows us to apply the Cox model to solving the explanation problems as a linear approximation of some unknown function of covariates by considering coefficients of the covariates as quantitative impacts on the prediction. As a result, Kovalev at al. [38] proposed an explanation method called SurvLIME, which deals with censored data and can be regarded as an extension of LIME 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 underlying SurvLIME, Utkin et al. [62] 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 [37] is called SurvLIME-KS. This model uses the well-known Kolmogorov-Smirnov bounds to ensure robustness of the explanation model to cases of a small amount of training data or outliers of survival data. This paper presents a method for producing counterfactual explanations for survival machine learning black-box models, which is based on the above results concerning with explanation of these models. The approach 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 survival functions or cumulative hazard functions. 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 the 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 till 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. The algorithm implementing the method is called as SurvCoFact. 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. 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. A heuristic global optimization method called Particle Swarm Optimization (PSO), proposed by Eberhart and Kennedy [34], can be applied to solving the counterfactual explanation problem. The method is a population-based stochastic optimization technique based on swarm and motivated by the intelligent collective behavior of some animals [67].
In summary, the following contributions are made in this paper: 1. A statement of the counterfactual explanation problem in the framework of survival analysis is proposed and a criterion for defining counterfactuals is introduced.
2. It is shown that the counterfactual explanation problem can be reduced to a standard convex optimization problem with linear constraints when the black-box model is the Cox model. This fact leads to an exact solution of the counterfactual explanation problem.
3. The PSO algorithm is applied to solving the counterfactual explanation problem in a general case when the black-box model may be arbitrary. 4. The proposed approaches are illustrated by means of numerical experiments with synthetic and real data, which show the efficiency and accuracy of the approaches.
The paper is organized as follows. Related work is 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 on the case of survival analysis. In Section 5, it is shown how to reduce the counterfactual explanation problem to the convex programming by the black-box Cox model. The PSO algorithm is briefly introduced in Section 6. 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. Concluding remarks can be found in Section 9.

Related work
Local explanation methods. Due to importance of the machine learning model explanation in many applications, a lot of methods have been proposed to explain black-box models locally. One of the first local explanation methods is the Local Interpretable Model-agnostic Explanations (LIME) [53], which uses simple and easily understandable linear models to approximate the predictions of black-box models locally. Following the original LIME [53], a lot of its modifications have been developed due to a nice simple idea underlying the method to construct a linear approximating model in a local area around a test example. These modifications are ALIME [57], NormLIME [2], Anchor LIME [54], LIME-Aleph [50], GraphLIME [31], SurvLIME [38], SurvLIME-KS [37]. Another explanation method, which is based on the linear approximation, is the SHAP [44,60], which takes a game-theoretic approach for optimizing a regression loss function based on Shapley values. A comprehensive analysis of LIME, including the study of its applicability to different data types, for example, text and image data, is provided by Garreau and Luxburg [22].
An important group of explanation methods is based on perturbation techniques [20,21,48,65], which are also used in LIME. The basic idea behind the perturbation techniques is that contribution of a feature can be determined by measuring how a prediction score changes when the feature is altered [17]. Perturbation techniques can be applied to a black-box model without any need to access the internal structure of the model. However, the corresponding methods are computationally complex when samples are of the high dimensionality.
Most explanation methods deal with the point-valued results produced by explainable black-box models, for example, with classes of examples. In contrast to these models, outcomes of survival models are functions, for example, SFs or CHFs. It follows that LIME should be extended on the case of models with functional outcomes, in particular, with survival models. An example of such extended models is SurvLIME [38].
Some counterfactual explanation methods use combinations with other methods like LIME and SHAP. For example, Ramon et al. [51] proposed the so-called LIME-C and SHAP-C methods as counterfactual extensions of LIME and SHAP, White and Garcez [71] 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 [40], 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 an interpretable counterfactuals in a lot of methods [15,16,42].
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 [23]. Hendricks et al. [27] proposed counterfactual explanations in natural language by generating counterfactual textual evidence. Counterfactual "feature-highlighting explanations" were proposed by Barocas et al. [6] by highlighting a set of features deemed most relevant and withholding others.
A lot of other approaches have been proposed in the last years, but there are no counterfactual explanations of predictions provided by the survival machine learning systems. Therefore, our aim of the presented work 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 a lot of real application problems, especially in medicine, reliability analysis, risk analysis. Therefore, the machine learning approach for solving the survival analysis problems allows improving the survival data processing. A comprehensive review of the machine learning models dealing with survival analysis problems was provided by Wang et al. [69]. One of the most powerful and popular methods for dealing with survival data is the Cox model [13]. A lot of its modifications have been developed in order 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 [61] presented a modification based on the Lasso method. Similar Lasso modifications, for example, the adaptive Lasso, were also proposed by several authors [36,73,76]. The next extension of the Cox model is a set of SVM modifications [35,72]. Various architectures of neural networks, starting from a simple network [18] proposed to relax the linear relationship assumption in the Cox model, have been developed [26,33,52,78] 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 [8] called the random survival forest (RSF) [32,45,68,74].
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, in order to explain predictions of these models.

Basic concepts
In survival analysis, an example (patient) i is represented by a triplet ( is the vector of the patient parameters (characteristics) or the vector of the example features; T i is time to event of the example. If the event of interest is observed, T i corresponds to the time between baseline time and the time of event happening, in this case δ i = 1, and we have an uncensored observation. If the example event is not observed, T i corresponds to the time between baseline time and end of the observation, and the event indicator is δ i = 0, and we have a censored observation. Suppose a training set D consists of n triplets (x i , δ i , T i ), i = 1, ..., n. The goal of survival analysis is to estimate the time to the event of interest T for a new example (patient) with feature vector denoted by x by using the training set D.
The survival and hazard functions are key concepts in survival analysis for describing the distribution of event times. The SF denoted by S(t|x) as a function of time t is the probability of surviving up to that time, i.e., S(t|x) = Pr{T > t|x}. The hazard function h(t|x) is the rate of event at time t given that no event occurred before time t, i.e., h(t|x) = f (t|x)/S(t|x), where f (t|x) is the density function of the event of interest. The hazard rate is defined as Another important concept is the CHF H(t|x), which is defined as the integral of the hazard function h(t|x), i.e., The SF can be expressed through the CHF as S(t|x) = exp (−H(t|x)).

The Cox model
Let us consider main concepts of the Cox proportional hazards model, [30]. According to the model, the hazard function at time t given predictor values x is defined as Here h 0 (t) is a baseline hazard function which does not depend on the vector x and the vector b; Ψ(x) is the covariate effect or the risk function; b = (b 1 , ..., b d ) is an unknown vector of regression coefficients or parameters. It can be seen from the above expression for the hazard function that the reparametrization Ψ( In the framework of the Cox model, the SF S(t|x, b) is computed as Here H 0 (t) is the cumulative baseline hazard function; S 0 (t) is the baseline SF. It is important to note that functions H 0 (t) and S 0 (t) do not depend on x and b.
The partial likelihood in this case is defined as follows: Here R j is the set of patients who are at risk at time t j . The term "at risk at time t" means patients who die at time t or later.

Counterfactual explanation for survival models: problem statement
We consider a definition of the counterfactual explanation proposed by Wachter et al. [66] and rewrite it in terms of the survival models.

Definition 1 ([66])
Assume a prediction function f is given. Computing a counterfactual z =(z 1 , ..., z d ) ∈ R d for a given input x =(x 1 , ..., x d ) ∈ R d is derived by solving the following optimization problem: where l(·, ·) denotes a loss function, which establishes a relationship between the explainable black-box model outputs; µ(·, ·) a penalty term for deviations of z from the original input x, which is expressed through a distance between z and x, for example, the Euclidean distance; C > 0 denotes the regularization strength.
The function l(f (z), f (x)) encourages the prediction of z to be different in accordance with a certain rule than the prediction of the original point x. The penalty term µ(z, x) minimizes the distance between z and x with the aim to find nearest counterfactuals to x. It can be defined as 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 [14] or by minimizing the distance between the counterfactual z and the class prototypes [42].
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 way is to separate groups of patients in accordance with their feature vectors. However, 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: Then the optimization problem (7) can be rewritten as follows: 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: The condition for "classes" of patients can be also written as Let us unite the above conditions by means of the function where the 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 belong to 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 Its minimization encourages to increase the difference f (z) − f (x). Taking into account (8) and (13), the entire loss function can be rewritten and the following optimization problem is formulated: In sum, the optimization problem (15) has to be solved in order to find counterfactuals z. It can be regarded as the constrained optimization problem: subject to ψ(z) ≤ 0.
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 solving. 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: 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) > 0), 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 .

The exact solution for the Cox model
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 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 (τ ) the SF S(τ |x) under condition of using the Cox model can be represented as follows: and Hence, the mean value is where µ j = t j+1 − t j > 0. Denote u = z T b and consider the function Compute the following limits: The derivative of π(u) is Note that s exp(u) 0,j exp(u) ≥ 0 for all j and u; µ j ln(s 0,j ) ≤ 0 for all j. Hence, there holds 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 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 the function ζ(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.
It follows from the above that a − < 0 ≤ b − and the function ζ(u) = 0 at a single point u − which can be numerically computed. 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 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, It can be solved by means of standard programming methods.

Particle Swarm Optimization
The PSO algorithm proposed by Eberhart and Kennedy [34] can be viewed as a stochastic optimization technique based on swarm. There are several survey papers devoted to the PSO algorithms, for example, [67,70]. 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 (u) where f : R n → R has to be minimized. The PSO is implemented in the form of the following algorithm:

Initialization (zero iteration):
• N particles u 0 k N k=1 and their velocities v 0 k N k=1 are generated; • the best position p 0 k = u 0 k of the particle u 0 k is fixed; • the best solution g 0 = arg min k f (p 0 k ) is fixed.
2. Iteration t (t = 1, ..., N iter ): • velocities are adjusted: where w, c 1 , c 2 are parameters, r 1 and r 2 are random variables from the uniform distribution in interval [0, 1]; • positions of particle are adjusted: • the best positions of particle are adjusted: • the best solution is adjusted: 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 as well as N iter should be as large as possible. Its upper bound depends only on the time that we can spend on iteration. We take N = 2000, N iter = 1000.
Particles are generated by means of the uniform distributions U with the following parameters: Velocities are similarly generated as:

Application of the PSO to the survival counterfactual explanation
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 > 0, δ 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 According to the PSO algorithm, the initial positions of particles are generated as So, the optimal solution can be found in the hyperparallelepiped X : 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 z closest,train = z ct = arg min j L(x j ).
Let us introduce a sphere B = B(x, R ct ) with center x and radius R ct = x − z ct 2 . The sphere can be partially located inside the hyperparallelepiped X or can be larger it. 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: Loop: over all features of z: The initial positions of particles are generated as follows: and the above restriction procedure is used for all points except the first one: u 0 k = RP (u 0 k ), k = 2, ..., N . Positions of particle will be adjusted by using the following expression: Initial values of velocities are taken as v 0 k = 0, k = 1, ..., N . Let us point out properties of the above approach: • The optimization results are always located in the set M.
• The "worst" optimal solution is z ct because the optimization algorithm remembers the point z ct at the zero iteration as optimal, and the next iterations never give the worse solution, if every initial position u 0 k starting from k = 2 is out of the feasible set of solutions, i.e., ψ(u 0 k ) > 0.

Numerical experiments
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. The outputs of the trained models in the testing phase are SFs.
2. In order 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 analyze the numerical results, the following schemes of their verification is proposed. When the Cox model is used as a black-box model, we can get exact solution (see Sec. 5). This implies that we can exactly compute the counterfactual z ver . Adding the condition that the solution belongs to the hyperparallelepiped 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 extremely complex 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 taken 10 6 . In fact, this approach can be regarded as the perturbation method with the exhaustive search.
When the black-box model is the RSF, then, in contrast to the black-box Cox model, the accuracy of the proposed method is defined by the relationship between two distances z ver − x 2 and z opt − x 2 , which will be denoted as A. The large values of A indicate that the PSO provides better results than the procedure of computing z ver by means of many generated points.

Initial parameters of numerical experiments with synthetic data
Random survival times to events are generated by using the Cox model estimates. An algorithm proposed by Bender et al. [7] for survival time data for the Cox model with Weibull distributed survival times is applied to generate the random times. The Weibull distribution for generation has the scale λ 0 = 10 −5 and shape v = 2 parameters. For experiments, we generate two types of data having the dimension 2 and 20, respectively. The two-dimensional feature vectors are used in order to graphically illustrate results of numerical experiments. The corresponding feature vectors x are uniformly generated from hypercubes [0, 1] 2 and [0, 1] 20 . Random survival times T j , j = 1, ..., N , are generated in accordance with [7] by using parameters λ 0 , v, b as follows: where ξ j is the j-th random variable uniformly distributed in interval [0, 1]; vectors of coefficients b are randomly selected from hypercubes [0, 1] 2 and [0, 1] 20 . The event indicator δ j is generated from the binomial distribution with probabilities Pr{δ j = 1} = 0.9, Pr{δ j = 0} = 0.1.
For testing, two points are randomly selected from the hyperparallelepiped X in accordance with two cases: θ = 1, condition (11), and θ = −1, condition (12). For each point, two tasks are solved: with parameter θ = 1 and parameter θ = −1. Parameter r is also selected randomly for every task.

The black-box Cox model
The first part of numerical experiments is performed with the black-box Cox model and aims to show how results obtained by means of the PSO approximate the verified results obtained as the solution of the convex optimization problem with objective function (16) and constraints (29). These results are illustrated in Figs. 1-2. The left picture in Fig. 1 shows how m(x) changes depending on values of two features x 1 and x 2 . It can be seen from the picture that m(x) takes values from 280 (the Similar results cannot be shown for the second type of synthetic data when feature vectors have the dimensionality 20. Therefore, we give them in Table 1 jointly with numerical results for the two-dimensional data. Parameters r ver and r opt in Table 1 are defined as respectively. They show the relationship between original margin r and margins r ver and r opt obtained by means of the proposed methods. In fact, values of r ver and r opt indicate how the obtained counterfactuals fulfil condition (11) or condition (12), i.e., conditions ψ(z ver ) = r − r ver ≤ 0, ψ(z opt ) = r − r opt ≤ 0.
The last three columns also display the relationship between z ver , z opt and x. In particular, the value of z ver − z opt 2 can be regarded as the accuracy measure of the obtained counterfactual.     3-4. In this cases, z ver is computed by generating many points (10 6 ) and computing m(z) for each point. The counterfactual z ver minimizes the distance z ver − x under condition ψ(z ver ) ≤ 0. It can be seen from the pictures that z ver is again very close to z opt .
Results of experiments with training data having two-and twenty-dimensional feature vectors are given Table 2. It can be seen from Table 2 that z opt is not worse (d = 2) or even better (d = 20) than z ver obtained by means of generating the large number of random points (see velues of A).

Numerical experiments with real data
We consider the following real datasets to study the proposed approach: Stanford2 and Myeloid. The datasets can be downloaded via R package "survival".
The dataset Stanford2 consists of survival data of patients on the waiting list for the Stanford In this dataset, we do not consider the feature "sex" because it cannot be changed. Moreover, we consider two cases for the feature "trt" (treatment arm), when it takes values "A" and "B". In other words, we divide all patients into two groups depending on the treatment arm. As a result, we have three datasets: Stanford2 and Myeloid-A and Myeloid-B.

The black-box Cox model
Since examples from the dataset Stanford2 have two features which can be changed (age x 1 and T5 mismatch score x 2 ), then results of numerical experiments for this dataset can be visualized, and they shown in Figs. 5-6. We again see that points z ver and z opt are close to each other. The same follows from Table 3 which is similar to Table 1, but contains results obtained for real data. If to consider values in the last column of Table 3 as the method accuracy values, then one can conclude that the method provides outperforming results. This means that the Cox model used as a black-box model accurately supports the dataset, and the PSO provides a good solution.
Results of experiments with the Cox model trained on datasets Stanford2, Myeloid-A and Myeloid-B are shown in Table 3. One can see that the proposed method also provides exact results for datasets Myeloid-A and Myeloid-B.

The black-box RSF
Results for the black-box RSF trained on the dataset Stanford2 are given in Figs. 7-8. We again see that z ver is close to z opt . Results of experiments with datasets Stanford2, Myeloid-A and Myeloid-B are shown in Table 4. We again see that values of A positive for all datasets. This means that the

Conclusion
One of the important difficulties of using the proposed method is to take into account categorical feature. The difficulty is that the optimization problem cannot handle categorical data and becomes a mixed integer convex optimization problem whose solving is a difficult task in general. Sharma et al. [58] proposed a genetic algorithm called CERTIFAI to partially cope with the problem and for computing counterfactuals. The same problem was studied by Russel [56]. Nevertheless, an efficient solver for this problem is a direction for further research. The same problem takes place in using the PSO. There are some modifications of the original PSO taking into account the categorical and integer features [11,39,59]. However, their application to the considered explanation problem is another direction for further research.
We have studied only one criterion for comparison of the SFs of the original example and the counterfactual. This criterion is the difference of mean values. In fact, this criterion implicitly defines different classes of examples. However, other criteria can be applied to the problem and to separating the classes, for example, difference between values of SFs at time moment. The study of other criteria is also an important direction for further research. Another interesting problem is when the feature vector is a picture, for example, a computer tomography image of an organ. In this case, we have a high-dimensional explanation problem whose efficient solution is also a direction for further research.