Assessing Fairness in the Presence of Missing Data
Year: 2021

Abstract:
Missing data are prevalent and present daunting challenges in real data analysis. While there is a growing body of literature on fairness in analysis of fully observed data, there has been little theoretical work on investigating fairness in analysis of incomplete data. In practice, a popular analytical approach for dealing with missing data is to use only the set of complete cases, i.e., observations with all features fully observed to train a prediction algorithm. However, depending on the missing data mechanism, the distribution of complete cases and the distribution of the complete data may be substantially different. When the goal is to develop a fair algorithm in the complete data domain where there are no missing values, an algorithm that is fair in the complete case domain may show disproportionate bias towards some marginalized groups in the complete data domain. To fill this significant gap, we study the problem of estimating fairness in the complete data domain for an arbitrary model evaluated merely using complete cases. We provide upper and lower bounds on the fairness estimation error and conduct numerical experiments to assess our theoretical results. Our work provides the first known theoretical results on fairness guarantee in analysis of incomplete data.

heading: 1 Introduction
text: Mounting evidence [20, 8, 26, 42, 43] has suggested that powerful machine learning algorithms can be unfair and lead to disproportionately unfavorable treatment decisions for marginalized groups. In recent years, there has been a growing body of research on addressing the unfairness and bias of machine learning algorithms [12].
Meanwhile, missing data are ubiquitous and present daunting challenges in real-world data analysis. Particularly, missing data, if not adequately handled, would lead to biased estimations and improper statistical inferences [29]. As such, analysis of incomplete data has been an active research area [4, 33]. More recently, there is also a growing recognition that missing data may have deleterious impact on algorithmic fairness. For example, in medicine, bias caused by missing values in electronic health records is identified as a significant factor contributing to unfairness of machine learning (ML) algorithms used in medicine that may exacerbate health care disparities [36, 20]. However, there has been little reported research on assessing fairness of statistical and machine learning models using datasets that contain missing values.
In the presence of missing data, one popular approach for analysis, particularly in biomedical studies, is to use only the set of complete cases, i.e., observations with all features observed, discarding incomplete cases. As a result, we can define two related yet distinct data domains, namely, the complete case domain and the complete data domain (see Section 2.1). In many biomedical applications,
∗This work is partly supported by NIH grant R01GM124111 and RF1AG063481. The content is the responsibility of the authors and does not necessarily represent the views of NIH.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
samples in the complete data domain are considered to be randomly drawn from the target population of interest; in other words, the complete data domain is defined by the target population, and the ultimate goal is to apply the trained models to the complete data domain. Under some missing data mechanisms (see Section 2.1), the complete case domain is biased for estimation of fairness in the complete data domain, i.e., for the target population. As such, a fair algorithm in the complete case domain may show disproportionate bias towards marginalized groups in the complete data domain.
The growing body of literature on algorithmic fairness has been primarily focused on two types of fairness definitions, group fairness and individual fairness [12]. Group fairness emphasizes that members from different groups (e.g. gender, race etc.) should be treated similarly, while individual fairness pays more attention to treatment similarity between any two similar individuals. In this work, we investigate group fairness in analysis of incomplete data. It has been noted that fairness definitions may not be compatible with one another in a sense that it is not possible to achieve fairness simultaneously under multiple definitions [19]. In binary classification problems, demographic (or statistical) parity [10, 17] is a fairness notion that has been mostly studied. It states that the predicted outcome should be independent of sensitive attributes. However demographic parity can cause severe harm to prediction performance when the response is dependent of sensitive attributes. As an alternative, disparate mistreatment [47] states that misclassification level (e.g., in terms of overall accuracy, false negative rate, false discovery rate) should be similar between two sensitive groups. Similarly, [24] proposed equalized odds, which requires both false positive rate (FPR) and false negative rate (FNR) to be the same between two groups. In the regression setting, fairness is usually associated with the parity of loss between two groups [1, 35, 16]. To fix ideas, we propose to use in this paper accuracy parity gap as the fairness notion in learning tasks including classification and regression. We consider the technique of re-weighting on assessing fairness of a given algorithm using complete cases, in which different complete cases are assigned different weights. We show that if the weights are properly chosen, such approach can mitigate the estimation bias induced by the difference of domains. It is worth noting that our results can be generalized to other fairness notions such as equal opportunity and prediction error parity with respect to mean square error.
Several existing works [18, 44, 48, 39, 30, 21] are related to our work, but there are a number of fundamental differences between these works and ours. [18, 44, 48] investigated the impact of missing data on the fairness of downstream prediction models. [39] investigated fairness across different domains and provided an upper bound of fairness in the target domain given fairness estimate in the source domain. But their work does not deal with missing data and associated challenges and does not consider the technique of re-weighting. In addition, they provided only upper bounds on transferring fairness. [30] empirically evaluated fairness in the presence of missing data. They did not observe consistent fairness results when comparing different methods for handling missing values such as imputation. We provide theoretical guarantees on assessing fairness of algorithms via analysis of incomplete data. [21] proposed a causal graph-based framework for modeling the data missingness to guide the design of fair algorithms. But they considered the case when the entire sample is missing, while we consider all three missing mechanisms.
Our contributions: This work offers four novel contributions. First, it provides new insights for algorithmic fairness in analysis of incomplete data, which has not yet been well-investigated from a theoretical perspective. Second, we characterize the role and impact of the missing data mechanism on correction for the data domain shift in fairness estimation through the analysis of incomplete data. Third, there has been limited theoretical work on domain adaptation problems when estimated weights are used to correct for domain shift. [13] considers the setting where true inverse probability weights are known, [41] considers the case when weights are estimated, but their analysis of algorithm’s generalization performance is compared with a clipped empirical risk and does not provide any lower bound. In this work we investigate the setting where the weights are estimated from the data based on correctly or incorrectly specified propensity score models under a given missing data mechanism. Besides, both upper and lower bounds for fairness estimation bias are provided. Fourth, while the estimand that is of primary interest in domain adaptation literature is the prediction accuracy, the estimand of our interest takes the form of absolute difference between two prediction accuracy terms. To obtain a tighter upper bound, we conduct proper relaxations and incorporate Bennett’s inequality in the proof of Theorem 1, which is more nuanced. In addition, novel proof techniques are also presented in the proof of Theorem 2 after applying Bernstein’s inequalities, which does not appear in existing works for establishing lower bounds in domain adaptation problems such as [13].

heading: 2 Problem Formulation

heading: 2.1 Preliminaries on missing data
text: The data structure and notation used in this paper mostly follow the convention in the missing data literature [29] and are summarized in Figure 1. Suppose we have a random sample of n observations from a target population of interest. If there were no missing values, each observation/case zi := {xi, yi} ∈ X × Y (i = 1, . . . , n) consists of predictors xi ∈ X and label (response variable) yi ∈ Y . Denote the complete data matrix by Z in which the ith row is denoted by zi. Some entries of Z are missing and we define the indicator for observing zij or not as rij = 1zij is observed, where zij is the j-th feature in zi. Denote the corresponding indicator matrix by R. Let z(1)i denote the components of zi that are observed for observation i, and z(0)i denote the components of zi that are missing for observation i. For example, consider the case when there are two predictors and one response; if only zi1 is observed, then z(1)i = zi1, z(0)i = (zi2, zi3). We then define the observed data Z(1) as the collection of the observed components from all n observations, { z(1)i, i = 1, . . . , n } and the
missing data Z(0) as the collection of all the missing components, { z(0)i, i = 1, . . . , n } .
There are three primary missing data mechanisms, namely, missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR) [29]. Data are said to be MCAR if the distribution of R is independent of Z. For MAR, the distribution of R depends on Z only through its observed components, i.e., R ⊥ Z(0)|Z(1). For MNAR, the distribution of R depends on the missing components of Z. We seek to investigate fairness guarantee under all three mechanisms.
In the presence of missing data, observation i is said to be a complete case if it is fully observed (i.e. z(1)i = zi). Let Ri := 1zi is fully observed denote the indicator of complete cases. We can then define two different data domains (distributions) and we use the two terms, domain and distribution, interchangeably in the remainder of the paper. The complete data domain, denoted by DT , is the distribution of zi in the joint space X × Y . The complete case domain, denoted by DS , is the distribution of observations that have all variables fully observed: zi|Ri = 1. There is an important connection between these two domains and the aforementioned Z(1) and Z(0). Z(1) contains the data in all the complete cases and the observed data in the incomplete cases. Combining Z(1) and Z(0) yields the complete data, which follows the distribution DT (Figure 1). Under MCAR, the distribution of the complete cases in DS is the same as the distribution of the complete data in DT . However under MAR or MNAR, the data distributions in DS and DT can be different. For example, if missingness depends on gender and other features associated with gender, then females may have a substantially higher proportion of missing values and the feature distribution in females in DS may be very different from that in DT . As a result, an algorithm that is trained to be fair in DS may not be fair when evaluated in DT , noting again that in practice we typically are more interested in fairness guarantee in DT . We can view this as a domain shift from DS to DT . This shift in the setting of incomplete data can be characterized by the true propensity score (a.k.a, the probability of being a complete case) [23]. To mitigate the fairness estimation bias caused by the shift, different weights can be assigned to the data in the fairness estimator. Missing data mechanism plays a vital role in obtaining the ‘right’ weight. For MCAR, uniform weights are equivalent to the
true weights. For MAR, one can train popular nonparametric models to obtain a consistent estimator for the true weights. For MNAR, it’s much more difficult to obtain consistent estimators for the true weights. Hence one of our contributions is to characterize the role and impact of the missing data mechanism on bias correction in fairness estimation through the analysis of incomplete data.

heading: 2.2 Fairness estimand
text: We consider learning tasks that use features x ∈ X to predict response y ∈ Y . Each observation also has a binary sensitive attribute A ∈ {0, 1}. We are interested in assessing fairness of a prediction model g : X → Y in the complete data domain DT . Let Ea(g) := ETa |g(x) − y(x)| denote the prediction error, where Ta represents that the expectation is taken with respect toDTa , the distribution of complete cases that belongs to sensitive group A = a (Figure 1). To derive our theoretical results, we consider the following fairness estimand.
Definition 1 (Accuracy Parity Gap). For a given prediction model g, the accuracy parity gap of g is ∆T (g) = |E0(g)− E1(g)|, where subscript T indicates that the fairness estimand is defined in the complete data domain DT .
This definition has close connections with various fairness notions proposed in the literature. In binary classification tasks where the response is binary: y ∈ {0, 1}, the notion accuracy parity has been used in [47, 19, 50], which requires that the prediction accuracy between two sensitive groups to be equal: P (g(x) = y|A = 0) = P (g(x) = y|A = 1). Accuracy parity gap in this case is the absolute value of difference between above two quantities |P (g(x) = y|A = 0)− P (g(x) = y|A = 1)|. For regression tasks where the response y takes continuous value, fairness constraints on loss difference between two groups are adopted in [16, 35, 1]. Accuracy parity gap under such setting can be regarded as the difference of the mean absolute error (MAE) loss between two sensitive groups.

heading: 2.3 Fairness estimator
text: To estimate ∆T (g), one can first estimate E0(g) and E1(g) using the set of complete cases in DS . However, the resulting estimator can be biased because of the difference between DT and DS under MAR and MNAR. To mitigate such estimation bias, one useful approach is to assign weight ω(zi, Ai) to observation zi with sensitive attribute Ai, and calculate the weighted sum over the complete cases to estimate E0(g) and E1(g). Specifically, we define the weighted empirical risk (prediction error) using the complete cases as Êa(g, ω) := 1∑n
i=1 I(Ai=a)Ri
∑n i=1 I(Ai = a)Riω(zi, Ai)|g(xi)− yi|,
where a ∈ {0, 1}. Here we assume there is at least one complete case observed for each sensitive group (i.e. ∑n i=1 I(Ai = a)Ri ≥ 1). Then the proposed fairness estimator is defined as follows.
Definition 2 (Fairness estimator from complete cases). Suppose the weights assigned to complete cases are given by ω. Then the fairness estimator for prediction model g in the complete data domain is ∆̂S(g, ω) = |Ê0(g, ω)− Ê1(g, ω)|, where subscript S indicates that the estimator is obtained from the data from the complete case domain DS .
An ideal choice of ω is the normalized inverse of the propensity score, which can effectively mitigate the bias caused by the difference of DT and DS . Mathematically, we let π(zi, Ai) := PT (Ri = 1|zi, Ai) denote the true propensity score (PS) model. In practice, we typically do not know the true propensity scores or the true distribution of complete cases, so we need to estimate the propensity scores and use the empirical distribution of complete cases. Various statistical and machine learning models can be used to estimate the propensity scores, such as logistic regression, random forest [7], support vector machines [15] and boosting algorithms [45, 32].

heading: 3 Main Results
text: In this section we provide the main theoretical results for the proposed fairness estimator ∆̂S(g, ω) in the complete data domain DT . Throughout, we assume that the weights are normalized in the complete case domain for both sensitive groups a ∈ {0, 1}: ESaω(z, a) = 1 and are bounded away from 0 and infinity. Here Sa represents that the expectation is taken with respect to DSa , the distribution of complete cases in the sensitive group A = a. Of note, if the weights ω(z, a) are independent of z, then ω(z, a) ≡ 1. It is commonly assumed in the domain adaptation literature that,
without loss of generality, y in the regression setting takes value inside interval [b1, b2] for some real numbers b1 and b2 [14, 49].
While our theoretical results are derived for APG in Definition 1, the framework of our proofs can be adapted to other fairness notions and the resulting forms of Ea(g). Take the example of a fairness notion for the regression problem, which is defined as the difference of Lp loss with 1 ≤ p <∞. The form of the covering number and the variance term VarSi(ω(z)|g(x)− y(x)|) in the analysis should be adjusted. For binary classification, our framework can be adapted to all the five measurements of disparate mistreatment (accuracy, false positive rate, false negative rate, false mission rate and false discover rate) mentioned in Section 2 of [47]. Notably, false negative rate parity is also known as Equal Opportunity proposed in [24].

heading: 3.1 An upper bound
text: Let B = supzω(z, A) < +∞ denote the upper bound of the weights, Dωa = ESaω(z, a)2 denote the second moment of weights. ωDSa denotes the distribution whose probability density function at z equals to ω(z, a)fSa(z), with fSa being the probability density function in DSa . We further let na denote the number of complete cases in group a ∈ {0, 1} defined by the sensitivity feature. Without loss of generality, throughout the paper we assume n0 ≤ n1, that is, sensitive group A = 0 is always the minority group. Theorem 1 below provides an upper bound of the fairness estimation bias.
Theorem 1. Assume that g is from a hypothesis classH with VC dimension d (pseudo dimension if in the regression setting) and y ∈ [0, 1]. Assume Dωa ≤ na/8 for both groups. Then, for any δ > 0, with probability at least 1− δ, the following inequality holds:∣∣∣∆T (g)− ∆̂S(g, ω)∣∣∣ ≤ ∑
a∈{0,1}
dTV (DTa ||ωDSa) + √√√√ BCd(na, Dωa , δ) na [ (1 + Dωa B ) log(1 + B Dωa )− 1 ] (1)
where dTV (DTa ||ωDSa) denote the total variation distance between DTa and ωDSa , and that
Cd(na, D ω a , δ) = log (d+1)(8e)d+1 δ + d 2 log na 2Dωa if g ∈ {0, 1} is a classification model log 4δ ( 8e d )d + 3d2 log na 3 √
2Dωa if g ∈ [0, 1] is a regression model and na ≥ d
Remark 1. If y ∈ [b1, b2], the upper bound in Theorem 1 would be multiplied by b2 − b1 under a new assumption Dωa (b2 − b1) ≤ na/8 for both sensitive groups. Remark 2. Dωa is always upper bounded by B2. In particular, for the MCAR mechanism, Dωa = 1.
The detailed proof of Theorem 1 is in Appendix A. Since the fairness estimand is defined based on prediction error, there is some similarity between the upper bound in Theorem 1 and the upper bounds for learning error in the domain adaptation [6, 37] and survival analysis [5]. Our upper bound is obtained by the triangle inequality and the detailed analysis of generalization error. In the first term of the upper bound, ωDSa is an approximation to the complete data domain in sensitive group a, and dTV(DTa ||ωDSa) can be viewed as the approximation error. It follows that a less accurate approximation of the complete data domain would lead to a looser upper bound on fairness estimation error. In the second term Cd(na, Dωa , δ) is proportional to log na, so the term is of order O(log n0/n0) 1 2 . It also follows that for a fixed total number of complete cases, the upper bound increases with sample imbalance between two groups defined by A. In addition, the missing data mechanism impacts the second moment of estimated weights Dωa ∈ [1, B2]. If a missing data mechanism leads to a larger second moment of the weights, the upper bound for the fairness estimation error would be looser. Furthermore, when the weights ω are defined using the true propensity scores, we call the resulting ω0(zi, Ai) = [π(zi, Ai)ES {1/π(z, A)|A = Ai}]−1 as the true weights. When true weights are adopted, the first term in the upper bound in Theorem 1 vanishes and we have the following result.
Corollary 1. If ω(zi, Ai) = ω0(zi, Ai), then ∆̂S(g, ω) is consistent for estimating ∆T (g).
There are several implications from Corollary 1. Under MCAR, setting ω(zi, Ai) = 1 would yield a consistent (unweighted) estimator. Under MAR and MNAR, since we typically do not know the true propensity scores π(zi, Ai), we replace π(zi, Ai) with its estimate π̂(zi, Ai) using a working model
which is subject to mis-specification. If a correctly-specified propensity score model is adopted, the second term in the upper bound would be the dominant term, in which the upper bound can be approximated from the observed data more easily. Otherwise, the first term would often be the dominant term, in which dTV(DTa ||ωDSa) can be even larger than that of the unweighted estimator.

heading: 3.2 A lower bound
text: We define σ2a(g, ω) := VarSa(ω|g − y|). Given ω is upper bounded by B and g, y ∈ [0, 1], we have that σ2a(g, ω) ≤ B2/4 by Popoviciu’s inequality on variance. We present the result on a lower bound for the fairness estimation error in Theorem 2. Theorem 2. If the weight ω(zi, Ai) is set to be ω0(zi, Ai) (true weights) and B2/σ2a(g, ω) ≤ n0, then the following hold with probability at least 71440 ,
12
√ σ20(g, ω0)
n0 + σ21(g, ω0) n1 ≥ ∣∣∣(E0(g)− E1(g))− (Ê0(g, ω)− Ê1(g, ω))∣∣∣ ≥ 1 24
√ σ20(g, ω0)
n0 + σ21(g, ω0)
n1 (2) Additionally, if ∆̂S(g, ω) ≥ 132 √ σ20(g,ω0) n0 + σ21(g,ω0) n1 , we have
∣∣∣∆T (g)− ∆̂S(g, ω)∣∣∣ ≥ 1 24
√ σ20(g, ω0)
n0 + σ21(g, ω0) n1
If ∆̂S(g, ω) ≤ 172 √ σ20(g,ω0) n0 + σ21(g,ω0) n1 , we have:
∣∣∣∆T (g)− ∆̂S(g, ω)∣∣∣ ≥ 1 72
√ σ20(g, ω0)
n0 + σ21(g, ω0) n1
The detailed proof of Theorem 2 is in Appendix B. The proof involves detailed analysis of the truncation probability for the fairness difference. Remark 3. If y is bounded in [b1, b2] in regression instead of the unit interval, Theorem 2 still holds. Remark 4. If the weight ω is based on the estimated π̂(zi, Ai) from a correctly specified propensity score model, then the optimal convergence rate for |ω(zi, Ai)− ω0(zi, Ai)| is Op((n0 + n1)−1/2). Based on the results from (1) and (2), the upper bound is of orderO((log n0/n0)1/2). With some additional assumptions on the data distribution, the lower bound can be shown to have orderO((n0)−1/2). (See Appendix B for a more detailed analysis) Remark 5. The bounds in Theorem 2 are established under fairly weak conditions on data distributions. If one is willing to make additional assumptions on the tail behavior (e.g., gaussian or sub-gaussian), the bounds would hold with higher probabilities. Remark 6. There is very limited work on the lower bound analysis in the fields of domain adaptation and fairness. The existing works [13, 27] have reported similar results as ours, i.e., the lower bound of generalization error holds with a low probability.
Our theoretical results have additional important implications related to the missing data mechanisms. First, as long as sample size is sufficiently large as required, Theorem 1 would always hold for all mechanisms. Under the MCAR mechanism, the true propensity score is a constant and hence can be regarded as known, and the results from Theorem 2 hold for the unweighted estimator. Under the MAR mechanism, the true propensity score is generally unknown. Under the MNAR mechanism, the propensity score model depends on missing values, so it cannot be estimated without making additional modeling assumptions. If the propensity score model is mis-specified under MAR or MNAR, the results in Theorem 2 are not applicable.

heading: 4 Numerical Experiments
text: In this section we empirically evaluate the bias of fairness estimation in both synthetic and real data sets, particularly the results in Theorem 1 and Theorem 2. As a reminder, our goal is to estimate fairness ∆T (g) (i.e., APG) defined in the complete data domain, while the estimator ∆̂S(g, ω) is obtained from the complete case domain.

heading: 4.1 Synthetic experiments
text: In our simulation experiments, we assess the upper bound in Theorem 1 in a classification task and the lower bound in Theorem 2 in a regression task. We further investigate the effects of different factors on the fairness estimation in regression tasks. In each experiment, we generate 10 predictors and a binary sensitive attribute A ∈ {0, 1} with n samples. Unless noted otherwise, the predictors are generated from Gaussian distributions: xij∼N (1− 2Ai, 0.52) (i = 1, . . . , n and j = 1, . . . , 10). We generate missing values among the last five predictors, from a pre-specified propensity score model. The propensity score model and responses y are set differently in different experiments. In all the experiments except the first one, we use a set of 2000 data (predictors are drawn from aforementioned distribution) to train a prediction algorithm g, where the dataset is balanced regarding the sensitive attribute. We use another set of data, which contains missing values, to calculate ∆̂S(g, ω). In our numerical experiments, it is difficult to obtain the analytical form of the true accuracy parity gap ∆T (g), so we approximate ∆T (g) using the Monte Carlo method with additional 100000 samples that are not used in obtaining ∆̂S(g, ω) or g.
Assessment of the upper bound in Theorem 1: In this experiment, we seek to answer the question: What practical guarantees can Theorem 1 provide for estimation of ∆T (g), the true fairness of a prediction model g in the complete data domain. Given ∆̂S(g, ω) and the right hand side of (1), we can construct an interval in which, with a high probability, ∆T (g) lies. We consider a binary classification problem with response yi ∼ Bernoulli((1 + exp(x>i β))−1), where β = (0.1, 0.1, 0.1, 0.1, 0.1, 1, 1, 1, 1, 1)>. Each sensitive group contains 50000 samples and the missingness of data is drawn according to the following propensity score model logit(π(zi, Ai)) = 0.25− 0.5Ai, where logit(p) = log p1−p . We use the complete cases among aforementioned samples to build a linear support vector machine g. Since the propensity score model only depends on sensitive attribute a, by definition, the true weights ω0(z, a) ≡ 1. We calculate ∆̂S(g, ω) from the complete cases with ω = 1 and calculate the right hand side of (1). It follows from Theorem 1 that with probability at least 95%, the true fairness ∆T (g) is covered by [0, 0.25]. We use additional 100000 data to approximate the APG value ∆T (g) ≈ 0.13. In this experiment setting, we observe that the upper bound in (1) is not significantly loose and might be able to provide useful information about the fairness of g in practice (e.g., g is not extremely unfair in the complete data domain).
Assessment of the lower bound in Theorem 2: In this experiment, we assess the upper and lower bounds on both regression and classification problems. Let β = (0.1, 0.1, 0.1, 0.1, 0.1, 1, 1, 1, 1, 1)> and = ( 1, . . . , n)> ∼ N (0, In). For regression problem, we consider yi = (x>i β)2 + i, for classification problem, we consider yi ∼ Bernoulli((1 + exp(x>i β))−1). We use linear SVM as the prediction model g, and the missingness of data is drawn from the following propensity score model logit(π(zi, Ai)) = −3 + 15 ∑5 j=1 xij . Of note, since Ri depends on only the fully observed features, the missing data mechanism is MAR. In this experiment, we use the true weights ω(z, A) = ω0(z, A) to calculate the fairness estimation bias |∆T (g)−∆̂S(g, ω)|. We check the assumptions and plot both
bounds in Figure 2. Observe that the lower bound from Theorem 2 is always smaller than the mean of fairness estimation bias, as the sample size n0 in the minority sensitive group changes from 102 to 104 while the ratio of n0/n is fixed at 1/2. Besides, in both problems, the lower bounds are close to the 5th percentiles of the fairness estimation bias, implying that they are not disproportionately loose. Our results from this and other unreported experiments lend support to that for certain data distributions such as those in this experiment, the lower bound may hold with a considerably higher probability than stated in Theorem 2. In addition, the slope of (average) fairness estimation bias is approximately the same as that of the lower bound, indicating the convergence rate of fairness estimation bias to be O(n−1/20 ) in this experiment setting. In addition, we notice that the confidence bands in Figure 2 seem to have equal length, given different sample sizes. However, since the y-axis has a log scale, the actual length of confidence band shrinks with increasing sample size.
Impact of different weights: We consider the same regression task in previous experiment, except that the predictors x are drawn from N (3− 6Ai, 0.52). Missingness of data is drawn under MAR using the following propensity score model logit(π(zi, Ai)) = − 15 ∑5 j=1 xij . We fix the sample ratio between two sensitive groups in the complete cases, which guarantees that both groups have the same number of complete cases (in expectation). In all the remaining synthetic experiments, we adopt random forest as the prediction model g. We calculate ∆̂S(g, ω) using seven different weights, namely, ω(zi, Ai) = 1 (i.e., unweighted estimator), the true weights ω(zi, Ai) = ω0(zi, Ai) and 5 weights obtained via ω(zi, Ai) = ω̂0(zi, Ai) := Riπ̂(zi,Ai) ∑n j=1{I(Aj=Ai)Rj/π̂(zj ,Aj)}
where π̂(zi, Ai) is estimated from two different logistic regression models (correctly and incorrectly specified), random forest (RF) [7], support vector machine (SVM) [15] and extreme gradient boosting (XGB) [11]. In particular we use {xij} to fit the first logistic regression model, which is correctly specified. In the second model, we use {x3ij}, which leads to an incorrectly specified logistic model. Figure 3-(a) shows that using the true weights or correctly specified logistic regression model yields similar performance and leads to smaller fairness estimation bias |∆T (g)− ∆̂S(g, ω)| than the incorrectlyspecified logistic model and the other propensity score models. At the first glance, it might seem surprising that XGB and random forest, as more expressive algorithms, do not outperform other propensity score estimators. A potential reason is that the true propensity score is the probability of observing all variables in a sample, which is often a number between 0 and 1 (bounded away from 0 by assumption). However, the data used to fit a propensity score model are binary labels (0 for missing data and 1 for observed data). The more expressive models such as XGB and random forest is more likely to predict label 1 for observed data and 0 for missing data accurately, but it does not always yield more accurate estimation of the true propensity score.
Impact of sample imbalance: We consider the same regression task defined in the previous experiment. We vary the total sample size n from 103 to 105, and for each fixed n we examine different levels of sample imbalance between the two sensitive groups by varying the ratio of n1/n0 from 1 to 20. We draw missingness of data under MAR from the following propensity score model, logit(π(zi, Ai)) = −1 + 15 ∑5 j=1 xij . To estimate the fairness ∆T (g), we use estimated propensity scores from the correctly specified logistic regression model. The resulting fairness estimation bias |∆T (g)− ∆̂S(g, ω)| is shown in Figure 3-(b). For a fixed n increasing sample imbalance leads to larger bias, suggesting that sample imbalance could harm the fairness estimation.
Impact of disparity in data distribution between sensitive groups: We consider the same regression task with n varying from 103 to 105, in which x is drawn from N (1 − 2MAi, 0.52) with M controlling the difference in data distribution between the two sensitive groups. Missingness of data is generated under MAR using the following model logit(π(zi, Ai)) = 2 − 4Ai. We use a correctly-specified logistic regression model to obtain the weights for estimating ∆̂S(g, ω). As shown in Figure 3-(c), the fairness estimation bias increases as M increases for a given n. This suggests that it can be harder to guarantee fairness in the complete data domain when the disparity in data distribution between the two sensitive groups become more pronounced.

heading: 4.2 Real data experiments
text: We conduct analyses of two real datasets, one from the COMPAS and the other from the ADNI. We consider the task of building prediction model and assess its fairness (in the target population) using the same dataset containing missing values. In each experiment, we randomly split the real dataset
into two subsets. In the first subset, we generate missing values, and the complete cases in this subset are used to train a random forest prediction model g and estimate its fairness in the complete data domain. The true fairness ∆T (g) is approximated using the entire second subset. Particularly, we generate missing values under three settings, namely, MCAR, MAR and MNAR respectively.
To assess the impact of weight specifications, we compare multiple options of ω, a) ω(zi, Ai) = 1 (i.e., unweighted), b) the true inverse probability weights ω(zi, Ai) = ω0(zi, Ai), and c) ω(zi) = ω̂0(zi, Ai) = Ri/(π̂(zi, Ai) ∑n j=1 {I(Aj = Ai)Rj/π̂(zj , Aj)}) where π̂(zi, Ai) is obtained from logistic regression, RF, SVM and XGB. To evaluate the impact of sample imbalance between two sensitive groups, we fix the total sample size and use logistic regression to compute weights ω̂0(zi, Ai). In all real data analyses, the logistic regression model is the correctly specified model for π(zi, Ai) under MAR mechanisms. Of note, while we can obtain fairness estimators (in the complete case domain) using real datasets with actual missing values, we would not be able to assess the bias in fairness estimation since the true fairness of an algorithm (in the complete data domain) would not be available. In addition, we would not be able to use the true propensity score model for missingness as a valuable benchmark. As such, we chose to generate artificial missing values in real datasets, which allows us to compute bias in estimating fairness and use the true propensity score model as a benchmark while making the experiments more realistic.
COMPAS recidivism dataset: Correctional Offender Management Profiling for Alternative Sanctions (COMPAS) [34] is a risk assessment instrument developed by Northpointe Inc. The dataset analyzed in this work contains records of defendants from Broward County from 2013 and 2014. Prior work has demonstrated the bias of predictions from COMPAS towards certain groups of defendants defined by race, gender and age etc. [2]. In our analysis, gender is treated as the sensitive attribute and nine numerical features are used to predict two-year recidivism (defined by arrest within 2 years) [38]. We generate missing values for the last feature and the outcome variable under three missing mechanisms: MCAR, logit(π(zi, Ai)) = 0.8; MAR, logit(π(zi, Ai)) = 3 + 2 ∑5 j=1 xij ; MNAR, logit(π(zi, Ai)) = −2y − 2xi9. As shown in Table 1, all options of ω lead to comparable results under MCAR, noting that all of them including the unweighted estimator are valid under MCAR.
Under MAR, the use of true weights leads to the least bias and followed by logistic regression, noting that logistic regression is the correctly specified model for π(zi, Ai) in this data analysis. Under MNAR, all working propensity score models are mis-specified and yield larger bias than the true weights. To study the impact of sample imbalance, we fix the total sample size as 800 and vary the proportion of samples in the two sensitive groups, nmale/nfemale, from 1 to 9. The curve of fairness estimation bias is shown in Figure 4. The results show that in more imbalanced data, fairness estimation bias is larger, consistent with our findings in the simulations.
ADNI gene expression data: The dataset from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) contains gene expression and clinical data for 649 patients. In our analysis, training set contains 500 samples and we only include the top 1000 transcriptomic features with highest positive correlation with gender, the sensitive feature. The outcome variable is the VBM right hippocampal volume. Missing values are generated for the last 900 features under the three missing data mechanisms: MCAR, logit(π(zi, Ai)) = 0.5; MAR, logit(π(zi, Ai)) = −2− 15 ∑10 j=1 xij ; MNAR,
logit(π(zi, Ai)) = −2− 15 ∑110 j=101 xij . As shown in Table 1, the main findings are consistent with those from the analysis of the COMPAS dataset. Regarding the effect of sample imbalance with a fixed total sample size 150, Figure 4 displays the curve of fairness estimation bias, showing the same patterns as in the analysis of the COMPAS dataset.

heading: 5 Discussions
text: This work provides the first known theoretical results on fairness guarantee in analysis of incomplete data. The bounds in Theorems 1 and 2 quantify the convergence rates of fairness estimation, which can help understand the impact of the sample size and guide choosing adequate sample size in practice. In addition, if the weights used in fairness estimation are not close to the true propensity score weights, the total variation distance term in the upper bound may become the dominant term. Thus, it is vital to use proper weights in fairness estimation. Our result can be adapted to other domain shift settings such as in the presence of selection bias where certain groups are under-represented due to biased sampling. A limitation of our work is that we only consider the analysis of complete cases. Another commonly-used alternative approach is to impute missing values first and then assess the fairness of the prediction model. However, there are several potential challenges for assessing fairness guarantee after imputation. Fairness estimation based on imputed data depends on the specific imputation method used, as different imputation methods have different operating characteristics and theoretical properties. In particular, the distributional properties of imputed data would play a vital role in understanding the theoretical properties of fairness estimation using imputed data. Such properties remain under-investigated for many popular imputation methods. For example, theoretical properties of multiple imputation via chain equations [9] and missForest [40] have not been well-established except for some very restrictive settings. Existing works on matrix completion (e.g. [31]) are focused on assessing imputation error, and there is little work on their statistical properties. Theoretical properties of popular deep learning imputation methods such as Mida [22], Gain [46] and misGan [28] have not been studied. We anticipate our work to lay a foundation for future research on assessing fairness in the presence of missing data.

Review summaries:
This paper discusses how to estimate fairness measures in the presence of missing values. Especially, it compares the fairness measure for the complete data domain (data distribution assuming no missing values) and the fairness measure estimated with reweighing on the complete case domain (data distribution after dropping all the data points that have missing values). They provide an upper bound and a lower bound on the difference of the two. They show experimental results of estimating the fairness measures on synthetic data and two real-world datasets to evaluate.  

In this paper the authors tackle the problem of ensuring fairness between subgroups in the presence of missing data, as the missingness can be biased between classes. The authors propose upper bound as well as lower bounds on the generalization error when propensity weights have to be estimated. The former is novel but classical, the latter is novel but also greatly original.

This paper provides theoretical guarantees (upper and lower bounds) on the estimation error of fairness, in the setting where missing data is present and only the complete rows are used for estimation. In particular, the authors consider the “accuracy parity gap” as the fairness metric, which can be applied to both classification and regression tasks. The paper also includes some numerical experiments to validate the theoretical findings.

The paper presents  theoretical results about fairness and missing data. The authors show theoretical bounds for the fairness metric "accuracy parity gap" considering weights. The weights are estimated from the data based on correctly or incorrectly specified propensity score models under a given missing data mechanism.  Some experiments on synthetic and real data (where some parameters are controlled) are useful to demonstrate the validity of  the included theoretical.