{'fraction_non_alphanumeric': 0.047476906942035414, 'fraction_numerical': 0.017990399741841797, 'mean_word_length': 4.162848812994586, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unit Selection with Nonbinary Treatment and Effect

The unit selection problem aims to identify a set of individuals who are most likely to exhibit a desired mode of behavior, for example, selecting individuals who would respond one way if encouraged and a different way if not encouraged. Using a combination of experimental and observational data, Li and Pearl derived tight bounds on the"benefit function", which is the payoff/cost associated with selecting an individual with given characteristics. This paper extends the benefit function to the general form such that the treatment and effect are not restricted to binary. We propose an algorithm to test the identifiability of the nonbinary benefit function and an algorithm to compute the bounds of the nonbinary benefit function using experimental and observational data.


Introduction
Several areas of industry, marketing, and health science face the unit selection dilemma. For example, in customer relationship management (Berson, Smith, and Thearling 1999;Lejeune 2001;Hung, Yen, and Wang 2006;Tsai and Lu 2009), it is useful to determine the customers who are going to leave but might reconsider if encouraged to stay. Due to the high expense of such initiatives, management is forced to limit inducement to customers who are most likely to exhibit the behavior of interest. As another example, companies are interested in identifying users who would click on an advertisement if and only if it is highlighted in online advertising (Yan et al. 2009;Bottou et al. 2013;Li et al. 2014;Sun et al. 2015). The challenge in identifying these users stems from the fact that the desired response pattern is not observed directly but rather is defined counterfactually in terms of what the individual would do under hypothetical unrealized conditions. For example, when we observe that a user has clicked on a highlighted advertisement, we do not know whether they would click on that same advertisement if it were not highlighted.
The binary benefit function for the unit selection problem was defined by Li and Pearl (Li and Pearl 2019) (we will call this Li-Pearl's model), and it properly captures the nature of the desired behavior. Using a combination of experimental and observational data, Li and Pearl derived tight bounds of the benefit function. The only assumption is that the treatment has no effect on the population-specific characteristics. Inspired by the idea of Mueller, Li, and Pearl (Mueller, Li, and Pearl 2021) and Dawid et al. (Dawid, Musio, and Murtas 2017) that the bounds of probabilities of causation could be narrowed using covariates information, Li and Pearl (Li and Pearl 2022c) narrowed the bounds of the benefit function using covariates information and their causal structure. However, the abovementioned studies are based on binary treatment and effect. Recently, researchers have shown interest in developing bounds for probabilities of causation with nonbinary treatment and effect. Zhang, Tian, and Bareinboim (Zhang, Tian, and Bareinboim 2022), as well as Li and Pearl (Li and Pearl 2022a), proposed nonlinear programming-based solutions to compute the bounds of nonbinary probabilities of causation numerically. Li and Pearl (Li and Pearl 2022b) provided the theoretical bounds of nonbinary probabilities of causation. The benefit function is a linear combination of probabilities of causation; therefore, in this paper, we focus on discovering the bounds of any benefit function without restricting them to binary treatment and effect.
Consider the following motivating scenario: a clinical study is conducted to test the effectiveness of a vaccine. The treatments include vaccinated and unvaccinated. The outcomes include uninfected, asymptomatic infected, and infected in a severe condition. The benefited individuals include the following: the individual who would be infected in a severe condition if unvaccinated and would be asymptomatic infected if vaccinated, the individual who would be infected in a severe condition if unvaccinated and would be uninfected if vaccinated, and the individual who would be asymptomatic infected if unvaccinated and would be uninfected if vaccinated. The harmed individuals include the following: the individual who would be asymptomatic infected if unvaccinated and would be infected in a severe condition if vaccinated, the individual who would be uninfected if unvaccinated and would be infected in a severe condition if vaccinated, and the individual who would be uninfected if unvaccinated and would be asymptomatic infected if vaccinated. All others are unaffected individuals. The researcher performing the clinical study has collected both experimental and observational data. The researcher then wants to know the expected difference between benefited and harmed individuals to emphasize the effectiveness of the vaccine.
We cannot apply Li-Pearl's model because we have two treatments and three outcomes. In this paper, we extend Li-Pearl's benefit function to general form without restricting them to binary treatment and effect. We will provide an algorithm to test the identifiability of the nonbinary benefit function and an algorithm to compute the bounds of the nonbinary benefit function using experimental and observational data.

Preliminaries
In this section, we review Li and Pearl's binary benefit function of the unit selection problem (Li and Pearl 2019), and the theoretical bounds of the probabilities of causation recently proposed by Li and Pearl (Li and Pearl 2022b).
In this paper we use the language of counterfactuals in structural model semantics, as given in (Galles and Pearl 1998;Halpern 2000). we use Y x = y to denote the counterfactual sentence "Variable Y would have the value y, had X been x". For simplicity purposes, in the rest of the paper, we use y x to denote the event Y x = y, y x to denote the event Y x = y, y x to denote the event Y x = y , and y x to denote the event Y x = y . we assume that experimental data will be summarized in the form of the causal efects such as P (y x ) and observational data will be summarized in the form of the joint probability function such as P (x, y). If not specified, the variable X stands for treatment and the variable Y stands for effect.
Individual behavior was classified into four response types: labeled complier, always-taker, never-taker, and defier. Suppose the benefit of selecting one individual in each category are β, γ, θ, δ respectively (i.e., the benefit vector is (β, γ, θ, δ)). Li and Pearl defined the objective function of the unit selection problem as the average benefit gained per individual. Suppose x and x are binary treatments, y and y are binary outcomes, and c are population-specific characteristics, the objective function (i.e., benefit function) is following (If the goal is to evaluate the average benefit gained per individual for a specific population c, argmax c can be dropped.): argmax c βP (y x , y x |c) + γP (y x , y x |c) + +θP (y x , y x |c) + δP (y x , y x |c).
Using a combination of experimental and observational data, Li and Pearl established the most general tight bounds on this benefit function (which we refer to as Li-Pearl's Theorem in the rest of the paper). The only constraint is that the population-specific characteristics are not a descendant of the treatment. Li and Pearl (Li and Pearl 2022b) provided eight theorems to compute bounds for any type of probabilities of causation with nonbinary treatment and effect. Suppose variable X has m values and Y has n values, the following probabilities of causation are bounded. Besides, if the probabilities of causation are conditioned on a population-specific variable c that is not affected by X, then all the theorems still hold (we provided the extended theorems from Li and Pearl in the appendix).
The benefit function is a linear combination of the probabilities of causation; therefore, we define the general benefit function for the unit selection problem based on Li and Pearl's results.

Counterfactual Formulation of the Unit Selection Problem
Based on Li and Pearl (Li and Pearl 2019), the objective is to find a set of characteristics c that maximizes the benefit associated with the resulting mixture of different response types of individuals. Let X denotes the treatment with m values and Y denotes the effect with n values. Therefore, we have n m different response types (i.e., one response type means assigning one effect to each of the treatments). Suppose the benefit of selecting an individual are (α 1 , ..., α n m ) (we call (α 1 , ..., α n m ) as benefit vector). Our objective, then, should be to find c that maximizes the following expression (If the goal is to evaluate the average benefit gained per individual for a specific population c, argmax c can be dropped): argmax c α 1 P (y 1x 1 , y 1x 2 , ..., y 1x m |c) + α 2 P (y 1x 1 , y 1x 2 , ..., y 2x m |c) + ... α n P (y 1x 1 , y 1x 2 , ..., y nx m |c) + ... α n m−1 +1 P (y 2x 1 , y 1x 2 , ..., y 1x m |c) + ... α n m P (y nx 1 , y nx 2 , ..., y nx m |c).
Note that c can be interpreted as the population-specific variable, the only assumption is that the treatment X has no effect on the population-specific variable. Recall from Li and Pearl's paper (Li and Pearl 2019), the benefit vector is provided by the decision-maker who uses the model.
In the next section, we will provide an algorithm that could check whether a given benefit function with the benefit vector is identifiable with purely experimental data (i.e., we can find the exact value of the benefit function rather than bounds). If it is not identifiable we will then provide an algorithm that computes the bounds of the benefit function given the benefit vector using experimental and observational data.

Identifiability of Benefit Function
Recall that in binary case, the conditions of identifiable are gain equality (i.e., β + δ = γ + θ) or monotonicity (i.e., P (y x , y x ) = 0) (Li and Pearl 2019). Here, it is complicated in nonbinary cases, therefore; we provide an algorithm to test whether a given benefit function with the benefit vector is identifiable with purely experimental data. Theorem 1. Suppose variables X has m values x 1 , ..., x m and Y has n values y 1 , ..., y n . Then the benefit function f (c) is identifiable if Algorithm 1 returns (True, res), and res is the value of the benefit function.
The correctness of the algorithm simply follow the fact that n m−1 terms P (..., y ix j , ...|c) = P (y ix j |c). Therefore, if there exist such n m−1 terms in the benefit function, then we can obtain an equivalent benefit function by replacing one of the n m−1 terms with experimental data P (y ix j |c). We exhausted all equivalent benefit functions to check if we could replace all the counterfactual terms with experimental data (i.e., identifiable).

Bounds of Benefit Function
If Algorithm 1 returns false, we then need to compute the bounds of the benefit function using experimental and observational data. We first obtain the bounds of the probabilities of causation, P (y 1x 1 , y 1x 2 , ..., y 1x m |c), ..., Algorithm 1: Check identifiability of the benefit function Input: a, the benefit function,where a[i] is a m + 1 tuple that stands for ith term in the benefit function. If the ith term is α i P (y i1 x1 , y i2 x2 , ..., y im xm |c), then e, the adjusted value of the benefit function.
where a[1, ..., n m ] corresponding to the original benefit function.
All lists in this algorithm start with index 1.
Output: (identifiable, value), a tuple, where identifiable = True if the given benefit function is identifiable and value is the value of the benefit function.
, then the sum of these n m−1 terms without coefficients is equal to P (y a[j1][r] xr ), we then recursively solve the equivalent benefit function. 14: for every n m−1 pair in a, 18:  Table 1: Experimental data of the clinical study. Here, 600 people were forced to take the vaccine and 600 people were forced to take no vaccine. P (y 1x 1 , y 1x 2 , ..., y nx m |c), ..., P (y 2x 1 , y 1x 2 , ..., y 1x m |c), ..., P (y nx 1 , y nx 2 , ..., y nx m |c), by Li and Pearl's theorems (Li and Pearl 2022b). We then have the following theorem.
Theorem 2. Suppose variables X has m values x 1 , ..., x m and Y has n values y 1 , ..., y n . Then the bounds of the benefit function f (c) is obtained by Algorithm 2.
Again, the correctness of the algorithm simply follow the fact that n m−1 terms P (..., y ix j , ...|c) = P (y ix j |c). We exhausted all equivalent benefit functions and take the maximum of all the lower bounds and take the minimum of all the upper bounds of equivalent benefit functions.

Example: Effectiveness of a Vaccine
Recall the motivating example at the beginning, a clinical study is conducted to test the effectiveness of a vaccine. The treatments include vaccinated and unvaccinated. The outcomes include uninfected, asymptomatic infected, and infected in a severe condition. The researcher of the clinical study has collected both experimental and observational data.

Task 1
The researcher wants to know the expected difference between benefited and harmed individuals to emphasize the effectiveness of the vaccine.
Let X denotes vaccination with x 1 being vaccinated and x 2 being unvaccinated and Y denotes outcome, where y 1 denotes uninfected, y 2 denotes asymptomatic infected, and y 3 denotes infected in a severe condition. The experimental and observational data of the clinical study are summarized in Tables 1 and 2. Based on the clinical study, the researcher of the vaccine claimed that the vaccine is effective in controlling the severe condition, the number of severe condition patients dropped from 213 to only 36.
Algorithm 2: Compute the bounds of the benefit function Input: a, the benefit function,where a[i] is a m + 1 tuple that stands for ith term in the benefit function. If the ith term is α i P (y i1 x1 , y i2 x2 , ..., y im xm |c), then lb, the lower bound of all possible terms obtained from Li-Pearl's theorems, where lb[(i 1 , i 2 , ..., i m )] is the lower bound of P (y i1 x1 , y i2 x2 , ..., y im xm |c).   Table 2: Observational data of the clinical study. Here, 1200 people were free to the vaccine. 953 people chose to take the vaccine and 247 people chose to take no vaccine. Now consider the expected difference between benefited and harmed individuals. Recall the benefited individuals include the individual who would be infected in a severe condition if unvaccinated and would be asymptomatic infected if vaccinated, the individual who would be infected in a severe condition if unvaccinated and would be uninfected if vaccinated, and the individual who would be asymptomatic infected if unvaccinated and would be uninfected if vaccinated. The harmed individuals include the individual who would be asymptomatic infected if unvaccinated and would be infected in a severe condition if vaccinated, the individual who would be uninfected if unvaccinated and would be infected in a severe condition if vaccinated, and the individual who would be uninfected if unvaccinated and would be asymptomatic infected if vaccinated. All others are unaffected individuals. In order to maximize the difference between benefited and harmed individuals; therefore, we assign 1 to benefited individuals, assign −1 to harmed individuals, and 0 to all others in the benefit vector. The objective function (i.e., benefit function) is then f (c) = 0P (y 1x 1 , y 1x 2 |c) + P (y 1x 1 , y 2x 2 |c) + P (y 1x 1 , y 3x 2 |c) − P (y 2x 1 , y 1x 2 |c) + 0P (y 2x 1 , y 2x 2 |c) + P (y 2x 1 , y 3x 2 |c) − −P (y 3x 1 , y 1x 2 |c) − P (y 3x 1 , y 2x 2 |c) + 0P (y 3x 1 , y 3x 2 |c).

Task 2
The researcher of the clinic study claimed that the individual who would be infected in a severe condition if unvaccinated and would be uninfected if vaccinated and the individual who would be uninfected if unvaccinated and would be infected in a severe condition if vaccinated should be twice important than other individuals. Based on the clinical study, the number of severe condition patients dropped from 213 to only 36; therefore, the vaccine should be effective for the virus. Now consider the expected difference between benefited and harmed individuals. The benefit vector should be the same except assigning 2 to the individual who would be infected in a severe condition if unvaccinated and would be uninfected if vaccinated and assigning −2 to the individual who would be uninfected if unvaccinated and would be infected in a severe condition if vaccinated.
We plug the estimates and the benefit function into Theorem 1, the Algorithm 1 returns true (i.e., identifiable by experimental data) with value −0.167. The benefit function can be simplified as follow: Thus, the expected difference between benefited and harmed individuals is exactly −0.167 per individual. We can conclude that the vaccine is still ineffective for the virus.

Simulated Results
In this section, we show the quality of the bounds of the benefit function obtained by Theorem 2 using four common benefit vectors.
First, we set m = 2 (i.e., X has two values) and n = 3 (i.e., Y has three values). We set the benefit vector to one of the most common ones, (0, 1, 1, −1, 0, 1, −1, −1, 0), which is to evaluate the expected difference between benefited and harmed individuals. We randomly generated 1000 populations where each population consists of different fractions of nine response types of individuals. For each population, we then generated sample distributions (observational data and experimental data) compatible with the fractions of response types (see the appendix for the generating algorithm). The advantage of this generating process is that we have the real benefit value (because we know the fractions of the response types) for comparison. Each sample population represents a different instantiate of the population-specific characteristics C in the model. The generating algorithm ensures that the experimental data and observational data satisfy the general relation (i.e., P (x, y|c) ≤ P (y x |c) ≤ P (x, y|c) + 1 − P (x|c)). For a sample population i, let [a i , b i ] be the bounds of the benefit function from the proposed theorem. We summarized the following criteria for each population as illustrated in  Besides, the average gap of the bounds, (bi−ai) 1000 , is 0.330, which is also small compared to the largest possible gap of 6.
Second, we set the benefit vector to another common one, (−1, 1, 1, −1, −1, 1, −1, −1, −1), which is to evaluate the expected difference between benefited and unbenefited (i.e., unaffected and harmed) individuals. We again randomly generated 1000 populations where each population consists of different fractions of nine response types. The data generating process and all other factors remain the same. We summarized the same criteria for each population as illustrated in Figure 2.
From Figure 2, it is clear that the proposed bounds obtained from Theorem 2 are a good estimation of the real benefit. The lower and upper bounds are closely around the real benefit and the midpoints are almost identified with the real benefit. Besides, the average gap of the bounds, (bi−ai) 1000 , is 0.6520, which is also small compared to the largest possible gap of 9.
Third, we set the benefit vector to another common one, (0, 1, 1, 0, 0, 1, 0, 0, 0), which is to evaluate the expected benefited individuals. We again randomly generated 1000 populations where each population consists of different fractions of nine response types. The data generating process and all other factors still remain the same. We summarized the same criteria for each population as illustrated in Figure 3.
From Figure 3, it is clear that the proposed bounds obtained from Theorem 2 are a good estimation of the real benefit. The lower and upper bounds are closely around the real benefit and the midpoints are almost identified with the real benefit. Besides, the average gap of the bounds, (bi−ai) 1000 , is 0.3284, which is also small compared to the largest possible gap of 3.
Lastly, we set the benefit vector to the last common one, (0, 0, 0, −1, 0, 0, −1, −1, 0), which is to evaluate the expected harmed individuals (we set the benefit vector to −1 because we want to minimize the harmed individuals). We again randomly generated 1000 populations where each population consists of different fractions of nine response types. The data generating process and all other factors still remain the same. We summarized the same criteria for each population as illustrated in Figure 4.
From Figure 4, it is clear that the proposed bounds obtained from Theorem 2 are a good estimation of the real benefit. The lower and upper bounds are closely around the real benefit and the midpoints are almost identified with the real benefit. Besides, the average gap of the bounds, (bi−ai) 1000 , is 0.3266, which is also small compared to the largest possible gap of 3.

Discussion
We have shown that the proposed theorems are a good estimation of the non-binary benefit function using examples and simulated studies. One may concern about the computation complexity of Algorithms 1 and 2. They are for sure in exponential time. However, the m and n (i.e., values of X and Y ) are usually small constant, therefore, we do not need to worry about too much.

Conclusion and Future Work
We demonstrated the formalization of the general benefit function with nonbinary treatment and effect. We provided the algorithm to compute the bounds of the general benefit function and the algorithm to check whether the benefit function is identifiable with purely experimental data. Examples and simulation results are provided to support the proposed theorems.
Future studies could assess the statistical properties of the proposed bounds. How tight would the bounds be? Does it sufficient to make decisions? Which data, experimental or observational, would affect the bounds more? How would the number of values in treatment and effect affect the quality of the bounds? Another future direction could be to improve the bounds using covariate information as Li and Pearl (Li and Pearl 2019) did for the binary benefit function.
Proof. Similarly to Theorem 1, lines 1 to 12 in Algorithm 2 simply compute the bounds of the given benefit function f (encoded as a in the algorithm). On line 24 of the algorithm, we recursively call on another benefit function f (encoded as na in the algorithm).

Li-Pearl's Bounds of Probabilities of Causation
The input, lb, ub, in Algorithm 2 depends on the bounds of probabilities of causation. The bounds of probabilities of causation recently proposed by Li and Pearl (Li and Pearl 2022b) is not conditional C. However, nothing is changed if conditioning on a variable C that is not affected by X. We listed the conditional version of the eight theorems proposed by Li and Pearl. The proof of eight theorems is exactly the same, except every probability should be conditioned on C.
Theorem 8. Suppose variable X has m values x 1 , ..., x m and Y has n values y 1 , ..., y n , and variable C is not affected by X, then the probability of causation P (y i1 xj 1 , ..., y i k xj k , x p |c), where 1 ≤ i 1 , ..., i k ≤ n, 1 ≤ j 1 , ..., j k , p ≤ m, j 1 = ... = j k = p, is bounded as following: , ..., y i k xj k |c)) +LB(P (y it xi t , x p |c)) − 1) , ..., y i k xj k |c)), where, LB(f ) denotes the lower bound of a function f and UB(f ) denotes the upper bound of a function f . The bounds of P (y i1 xj 1 , ..., y it−1 xj t−1 , y it+1 xj t+1 , ..., y i k xj k |c) are given by Theorem 7 or experimental data if k = 2 and the bounds of P (y it xi t , x p |c) are given by Theorem 5.

Distribution Generating Algorithm
Here, the sample distribution generating algorithm in the simulated studies is presented. It generated both experimental and observational data compatible with the fractions of response types of individuals. The data satisfy the general relation between experimental and observational data. Note that all four simulated studies shared the same distribution generating algorithm but with different benefit vectors. //rand(0, 1) is the function that random uniformly generate a number from 0 to 1. //f is the fractions of response types of individuals, f [0] = P (y 1x 1 , y 1x 2 |c k ), ..., f [8] = P (y 3x 1 , y 3x 2 |c k ). P (x 1 , y 1 |c k ) = rand(0, P (y 1x 1 |c k ));

33:
P (x 2 , y 3 |c k ) = P (x 2 |c k ) − P (x 2 , y 1 |c k ) − P (x 2 , y 2 |c k ); 34: //Validate the data, the experimental data and observational data should satisfies the following: P (x, y|c k ) ≤ P (y x |c k ) ≤ P (x, y|c k ) + 1 − P (x|c k ). if P (y ix j |c k ) < P (x j , y i |c k ) or P (y ix j |c k ) > P (x j , y i |c k ) + 1 − P (x j |c k ) then