Robust Propensity Score Computation Method based
on Machine Learning with Label-corrupted Data
Chen Wanga , Suzhen Wangb,1 , Fuyan Shib , Zaixiang Wangb

arXiv:1801.03132v1 [stat.ME] 9 Jan 2018

a

Dept. of Computer Science, University College London, Gower Street, London,
England, United Kingdom WC1E 6BT
b
Dept. of Health Statistics, Weifang Medical College, No.7166 BaoTong West Street,
Weifang, Shandong, China 261053

Abstract
In biostatistics, propensity score is a common approach to analyze the imbalance of covariate and process confounding covariates to eliminate differences
between groups. While there are an abundant amount of methods to compute propensity score, a common issue of them is the corrupted labels in
the dataset. For example, the data collected from the patients could contain samples that are treated mistakenly, and the computing methods could
incorporate them as a misleading information. In this paper, we propose
a Machine Learning-based method to handle the problem. Specifically, we
utilize the fact that the majority of sample should be labeled with the correct instance and design an approach to first cluster the data with spectral
clustering and then sample a new dataset with a distribution processed from
the clustering results. The propensity score is computed by Xgboost, and a
mathematical justification of our method is provided in this paper. The experimental results illustrate that xgboost propensity scores computing with
the data processed by our method could outperform the same method with
original data, and the advantages of our method increases as we add some artificial corruptions to the dataset. Meanwhile, the implementation of xgboost
to compute propensity score for multiple treatments is also a pioneering work
in the area.
Keywords: Biostatistics, Propensity Score, Spectral Clustering, Data
∗

This research is granted by China National Natural Science Foundation, Grant:
81473071
1
Corresponding Author

Sampling, Xgboost, Machine Learning
1. Introduction
Confounding covariates, or alternatively named as noise features, is a significant problem in studying Biostatistics data. In practice, data in this field
is usually collected with detailed information, thus it will contain some irrelevant and redundant features. For example, in a dataset of a certain disease,
the information about the marital status might be irrelevant to the cause of
the disease, and it will be acting as a confounding (noise) covariate. The
confounding covariates could bring negative effects in analyzing the data: it
could provide misleading information to reveal a false correlation between
variables (features) and labels, and it could introduce unnecessary statistical
differences between groups and lead to a incomparable situation. Consequently, recognizing and balancing confounding covariates have become an
open task in biostatistics research.
A common approach to address this problem is propensity score [1]. The
key idea of propensity score is to compute the ’likelihood’(propensity) of a
sample to be belonging into a certain group (usually termed as ’treatment’ in
biostatistics). And by matching the samples in different groups with similar
propensity scores and/or weighting the samples with propensity score, we
could reduce or eliminate the bias of confounding features. [2] experimentally shows that for some randomly generated data, propensity score could
effectively eliminate the statistical significance between groups. There are
multiple methods to compute propensity score, ranging from simple logistic
regression to cutting-edge Machine Learning methods. From the perspective
of Machine Learning, the task of propensity score computing is essentially
to train and utilize a classifier with probability as retrievable outputs (like
soft-max output). Consequently, although there are a variety of Machine
Learning algorithms which could be plunged to the task, the basic procedures are almost the same.
While the state-of-the-art Machine Learning methods are promising in this
area, they typically suffer from one problem: the corrupted data. Data is
typically collected from medical practice, and there could be some mistreated
samples in each group. The conventional approach to solving this problem
would be to hire some experts to select the correctly-labeled samples and
use them as a subset to train the classifier. However, this procedure will be
2

expensive in both time and money.
To tackle the issue with a smarter approach, we propose a novel Machine
Learning-based propensity score computation procedure in this paper. We
exploit the fact that the majority of samples should be treated (classified)
with the correct group (label) and patients (samples) be treated with different groups should have distinct data manifolds. Guided by this premise, we
design a system that first clusters the data with spectral clustering method,
and then computes the weight of the data by computing the proportion of
the number of samples through different classes (treatments). The weights
are processed with an interpolation method to perform a ’moderate softmax’
computation and it will lead to a distribution over clusters after being normalized. Based on this new distribution, we will sample a subset for the
usage of the training procedure. And finally, the propensity score is computed based on Xgboost proposed by Chen and Guestrin [3].
The rest of the paper is arranged in the following order: the second section
will briefly review literatures related to our research; the third part will introduce the proposed system and explain the Machine Learning algorithms
applied; the fourth portion will be mathematically analyzing the assumptions and properties of the proposed method; the experimental result for the
proposed method based on SEER data is illustrated in the fifth section, and
the outcomes are analyzed and compared; and finally, the last section of the
paper concludes our research and discusses future research for this topic.
2. Related Work
Despite the solid theoretical foundation of conventional Logistic Regression method in computing propensity score, efforts to utilize novel Machine
Learning algorithms in computing propensity score has emerged for a long
time. Huge research potentials of propensity score could be found in the
state-of-the-art algorithms, which could significantly outperform the traditional Logistic Regression. Existed papers related to this topic mainly focus
on implementing Machine Learning methods with advanced models (such
as SVM and Random Forest) and comparing different approaches with certain datasets. [4] tested different Machine Learning methods in computing
propensity score and demonstrated that boosting methods and tree algorithms could outperform Logistic Regression in general. [5] studied the theoretical foundations of implementing Machine Learning methods for dichotomous and multicategory outcome in biostatistics, and analyzed algorithms
3

like kNN, Random Forest and SVM. [6] compared the propensity score results
of Classification Tree Analysis (CTA) and Logistic Regression and argued
that although Logistic Regression method still has the lowest average standardized difference, CTA could provide the greatest predictive accuracy thus
it has strong potentials and could provide an alternative way in propensity
score computation. And for the tree-based method specifically, [7] compared
different classification and regression and discussed their performances.
There are also published literatures concentrating on improving the model
complexity to avoid potential model misspecified problem. A typical example
of these algorithms is the so-called ’super-leaner’ algorithm, which is roughly
equivalent to stacking algorithms in Machine Learning. The method was
originally proposed by [8] to process data with complex relationship between
input and output through a weighed combination of different learners. [9]
shows experimentally that the super learner method could achieve a significantly better performance under the situation of severe model misspecification. [10] proposed an algorithm to use a convex combination of generalized
propensity score computed by Super Learner. The authors then takes traumatic brain injury as a study case and discussed the statistical significance
between time for the patient to be transferred from emergency to specialist
and the treatment effect. [11] further studied the Super Learner method in
the application of electronic healthcare data.
Whilst there are a large amount of work discussing advanced algorithms in
computing propensity score, approaches dealing with propensity score under
corrupted data is relatively under-developed. [12] analyzed the problem of
propensity score with missing data and suggested to use multiple imputation
method. Similar with this paper, [13] analyzed propensity score computation
under missing covariate with multiple imputation and specifically solved the
issue with general location mixture models. To the best our knowledge, hitherto, there has not been any significant publication in studying propensity
score with mislabeled data, although this is common in bioscience. In the
area of Machine Learning, this issue could be regarded as corrupted label
problem, and there are some related literatures building algorithms for this
task [14][15]. However, in this work, we do not follow the idea of advanced
algorithmic target functions employed in these methods. Instead, we used
a clustering and normalization based method to tackle the problem. Based
on our assumption, the method should be valid for biostatistic data and the
rationality of the method is shown mathematically in our paper.

4

3. Methodology
3.1. Problem Setup
Assuming we have data (X, Y ) that generated from space X × Y, where
X is the input space and Y is a discrete output space. We now have m data
pairs (X1 , y1 ), (X2 , y2 ), ..., (Xm , ym ) drawn from a certain distribution. Then,
there will be an underlying conditional distribution:
p(y|X; θ)

(1)

And to infer the probability of each data point, we need to obtain parameter
θ that determines the model. With proper hypothesis space, the conventional
MLE of θ should be:
∗

θ = argmax
θ

m
Y

p(yi |Xi ; θ)

(2)

i=1

However, here our situation is the set Y has been replaced by Ŷ , which
contains many noisy(incorrect) labels. Then directly learning from (X1 , yˆ1 ),
(X2 , yˆ2 ), ..., (Xm , yˆm ) will incorporate some noise. Thus, it demands an approach to solve the problem.
Given the research context of medical and biological science, we could assume
the probability for a sample to be labeled (treated) correctly is higher than
the reverse situation. Mathematically, this assumption could be noted as:
p(ŷ = a|y = a) > p(ŷ 6= a|y = a), ∀a ∈ Y

(3)

And under this assumption, we could try to sample a subset of (X1 , y˜1 ),
(X2 , y˜2 ), ... , (Xk , y˜k ), k 6 m to fit the model, and find the parameter with:
θ̃ = argmax
θ

k
Y

p(ỹi |Xi ; θ)

(4)

i=1

that could minimize the loss L(θ̃, θ∗ ) between the estimated parameter and
the parameter estimated under uncorrupted data. However, since the true
data set is not visible for us, it is impossible to directly optimize over L(θ̃, θ∗ )
function. An alternative way would be to sample from a distribution that
could produce a subset close to the true dataset, and this is the core idea
of our noise-eliminating procedure. The next subsection will be introducing
this method and the employed Machine Learning algorithm.
5

3.2. Avoid Noisy Label: Sampling with Spectral Clustering and Normalization
Now we need a method to sample a subset of data with the distribution
as close as possible to the original correct samples. The method should have
the following properties:
• The proportion of mis-specified labels should monotonously decrease
as the performance of the algorithm increase
• The result should not be affected by the proportion of data in each
class (treatment)
• The results of the algorithm should be interpretable, and the interpretation should have a fixed methodology
To satisfy the above principles, we could assume some special properties for
the data. Without losing generality, we could assume the data, which is
classified (treated) with d different classes, comes from n underlying patterns under condition n > d, and each pattern will have a proportion of
pi , i = 1, 2, ..., n to be classified into class j ∈ {1, 2, .., d}. Formally, we could
specify the number of true samples under each class j:
Ntrue (j) =

n
X

pi N (i)

(5)

i=1

Where N (·) denotes the number of samples. In practice, we choose the number of n as n = d for a better interpretability, and this setting will be further
discussed in section 5 with the example of SEER data. Under this assumption, we will be able to sample the data with their underlying patterns, which
could be represented by clusters under the skeleton of unsupervised learning. Intuitively, an effective method would be to sample from the clustering
results, and now the question is which kind of technique should we use for
clustering.
The commonly-used clustering method is K-means and its variations. This
branch of algorithms are straightforward to implement and the computational complexities of them are generally acceptable for ordinary computers.
However, K-means related algorithms suffer from several problems: it is sensitive with initial clustering centers, and in practice we usually need to run
6

K-mean for several times with different initial clustering centroids to obtain
the best performance; it is vulnerable towards outliers, as few outliers will
make significant effects on the result; and most importantly, the clustering
results of K-means branch algorithms have a strictly spherical curvature, thus
it could only deal with simple data manifold, but in bio-medical science the
data manifold could be of great complexity.
A better alternative would be Spectral Clustering. Instead of measure the
Euclidean distance between samples, which will inevitably loss some manifold information, Spectral Clustering method utilize the Graph Distance that
could preserve local neighborhood informations. In Spectral Clustering, the
affine matrix that denotes the Graph Distance and its corresponding Graph
Laplacian will be computed, and then eigenvalue decomposition will be performed to capture the information in the so-called embedded space and the
clustering result could be obtained by applying ’ordinary’ clustering methods on this processed matrix. Spectral Clustering works especially well when
the number of clusters is not large, and it could deal with complicated data
geometries. The core problem of Spectral Clustering is the Neighbor Graph
(Similarity Graph), which represents the ’similarity’ between different data
points. Common methods to be implemented include -Neighbourhood, knearest Neighbourhood and fully-connected graph with positive semi-definite
kernels. One popular similarity graph matrix construction method is Gaussian Kernel method, which computes the similarity between data points xi
and xj with:
||xi − xj ||22
)
(6)
2σ 2
Which σ is a parameter that could be specified or computed by the following
equation:
q
Kg (xi , xj ) = exp(−

σ=

T

median(I(||X̂ − X̂ ||2 ))

(7)

Where X̂ is a m × m matrix with X̂ ij = ||xi ||22 , ∀j ∈ 1, 2, ..., m (duplicated
column matrix with each row as the square l2 norm of the data point) and I
stands for the operation to remove the diagonal values of the matrix.
The method could be computed with high efficiency for high dimension data.
However, it could be of great instability if some of the features are discrete
values. To improve the clustering method specifically for our data with a
mixture of continuous and discrete features, we proposed a feature-wise Similarity Matrix method here. For each pair of samples (xi , xj ), we compute
7

Data: Input data, separated by continuous features and discrete
features
Result: Clustering result for each data point
Compute the affine matrix of continuous features based on equation 6 ;
for Each Discrete Feature do
Compute the affine matrix of discrete features based on equation 9
and sum them with equation 8 ;
end
Perform Spectral Clustering with the affine matrix ;
Return the clustering label for each data ;
Algorithm 1: Spectral Clustering of the Data
the sum of the similarity across each feature, which could be denoted by:
K(xi , xj ) =

P
X

Kd (xip , xjp ) + Kg (x∗i , x∗j )

(8)

p=1

Where P is the number of discrete features(co-variates) and Kg (·, ·) stands
for Gaussian Kernel for continuous features computed by equation 6. For the
discrete features, we compute each of the similarity matrix with Delta Kernel:
(
1
, x i = xj
(9)
Kd (xi , xj ) = N (x)
0,
otherwise
Where N (x) is the number of samples of the specific class. And since all the
components of the similarity is positive semi-definite, the resulting matrix
would also be positive semi-definite. The Process of the detailed algorithm
could be shown in algorithm 1.
The practical implement of this method will be discussed in section 5. And
after clustering, we could further divide each cluster with sample in different
classes(treatments) and get a d × n matrix:


w11 w12 · · · w1n
w21 w22 · · · w2n 

W =
(10)
···

wd1 wd2 · · · wdn
Where d is the number of classes and n is the number of underlying patterns
(clusters).
8

Intuitively,now we would be able to choose the top k clusters at each class(treatment)
and use them as the new training set. However, the number of samples in
each cluster, which serves are an indication of underlying patterns, could be
imbalance and misguiding in our circumstance. A smarter way would be
to normalize the quantity of each cluster across the classes, which could be
representation as:
wij
∗
wij
= Pd
i=1 wij

(11)

This will give us the exact information about ’how large is the portion of
this pattern being assign with treatment i’, and it could best reflect and
satisfy the assumption denoted by equation 5. And here we need to augment the probability of the cluster with largest probability because we have
mis-specified data in each treatment class. Based on the assumption of 3,
the augmented probabilities should be belonging to the clusters that should
contribute to the treatment. And inspired by [16], we use an interpolation
∗
to get ŵ:
method to combine the first and second order value of wij
∗
ŵ = wij
+

γ
∗2
wij
1−η∗

(12)

Parameters γ denotes the trade-off between first and second order interpolation, and  is the proportion of corrupted labels (default as 0 even with some
mild corruptions, and we assume  ≤ 0.5) and η is a parameter denoting
how will we treat the importance of data corruption. Here we also have a
constraint to η for 0 ≤ η ≤ 2. In our program we set γ = 0.7 and η = 2. And
after the normalization, we could sample a new dataset with the proportion
of ŵ in each row of the matrix. Formally, it could be denoted by p∗ij :
wˆij
p∗ij = Pn
j=1

wˆij

(13)

Equation 13 is the equation we finally used in computing the probability
to sample a subset from each treatment class. In section 4 we will discuss
the rationality of it mathematically, but in this section, we will continue on
discussing the technique to compute propensity score.

9

3.3. Propensity Score Computing with Xgboost
As it is stated in section 2, implementing Machine Learning techniques
in computing Propensity Score has been studied for a long time. Propensity Score computing is very close to the probabilistic output in supervised
classification problem, thus a large amount of logistic or softmax functionbased approach could be introduced into the Propensity score area. Previous
studies generally suggest that boosting method is one of the most promising
branches of ML algorithms in computing propensity score. And in our work,
we use Xgboost as a novel Machine Learning method (designed in 2016) to
compute propensity score.
Xgboost is actually not a newly-proposed algorithm. Instead, it is built
on the algorithm of Gradient Boost Machine (GBM), and the novelty of
this platform is that it considered Structural Learning scores and designed a
better interface for Machine Learning considerations. The platform has R,
C++, Python, Scalar and Java API, and it could be conveniently modified
by simply change some parameters and arguments. It support softmax output for multiple-class probability, but using softmax could be too extreme.
The package also provide another mode called ’softprob’ model, which could
output the normalized probabilities of each class. And this is also the method
we adopted in this work.
Key parameters of Xgboost also include l2 regularization parameter, named
’lambda’, l1 regularization coefficient, named ’alpha’, and tree construction
and growth method, which is often selected automatically by heuristic search.
Meanwhile, since we are asking the Xgboost to output multi-class result, here
we use multi-class log loss to denote our training and validation loss. The
loss could be expressed by:
l(Y , P ) =

m
X

y Ti log(pi )

(14)

i=1

Where y i is the n × 1 one-hot vector that has value 1 at the true output
position and 0 otherwise; pi is the probability output vector. Y and P are
the m × n matrices that denotes all the data.
Xgboost is a well-developed platform which is ’efficient, flexible and portable’
(from their Github page). Here, we will not elaborate on the basic principles
of the GBM algorithm it utilizes. One who interested in could refer to the
original publication of the method.

10

3.4. Data Manifold Analysis with t-SNE
Another point we might interest in is the visualization of the data, which
could bring significant benefits for result interpretation. Visualizing the projection of the Manifold of data could provide us information about how data
in different classes could be overlapped or separated. And sometimes it could
even help us determine whether there exists significant difference between different groups (If there is not the experiment is probably unnecessary).
Typical data visualization approaches often project the data manifold to
a 2-d plane, which could be straightforwardly analyzed by human. Popular
methods used to compute data visualization include PCA(and pPCA, kPCA,
etc.), Isomap, Sammon Mapping, SNE(tSNE) and auto-encoder. The amjor
drawbacks of PCA and mapping methods are that they could not preserve
local structures and will thus lose the information of the data manifold in
the higher dimension. On the contrary, SNE(tSNE) could manipulate the
graph neibourhood structure. It is similar with Spectral Clustering method
we mentioned before, the difference is that here we compute the data point
in the lower dimension instead of clustering the Graph Laplacians.
To begin with, we will start to introduce SNE method. SNE stands Stochastic Neighbour Embedding, and it start with the transition matrix between
different data points in the original space, computed by:
(
exp(−||(xi )−(xj )||22 /(2σ 2 ))
P
i 6= j
2
2
j=1,j6=i exp(−||(xi )−(xj )||2 /(2σ ))
(15)
pi|j =
0
i=j
Similarly, we could compute the transition probabilities in the lower-dimension
space, which could define:
(
exp(−||(x∗i )−(x∗j )||22 /(2σ 2 ))
P
i 6= j
∗
∗
2
2
j=1,j6=i exp(−||(xi )−(xj )−||2 /(2σ ))
qj|i =
(16)
0
i=j
And now our purpose is to fit the model to get the data x∗ that could preserve the structural information. SNE is a method that consider to minimize
the KL divergence between pi|· and q·|i , and the target function could be describe as:

11

l((x∗ ), (x)) =

X

KL(pi|· |q·|i )

i

=

XX
i

pi|j log(

j

pi|j
)
qj|i

(17)

Equation 17 could be solved by gradient descent. And t-SNE is a variate
of SNE to modify the ’unsymmetrical’ KL-divergence. It defines a ’joint’
distribution of pi,j with:
pi|j + pj|i
(18)
2m
Where m is the number of samples. The similar method could be used on
q, and then we would have a ’symmetric’ target function. In this work, we
utilized the tSNE package in sk-learn API and specified certain parameters
to compute an optimal solution. Notice that in our situation features are
partially discrete and partially continuous, thus we compute the distance between samples explicitly with Euclidean for continuous features and matching
numbers for discrete features as the following equation:
pi,j =

N − Nm
(19)
N
Where N is the number of features(covariates) Nm is the number of features
that match in sample xi and xj .
Computing continuous and discrete features separately could cause another
problem: the imbalance of distance. Since the value of Dd (xi , xj ) has a maximum value of 1, where the distance of continuous data is not constrained,
the effect of Dd (xi , xj ) could be under-valued. An approach to deal with
this problem is to multiply a correction coefficient τ to Dd with value:
Dd (xi , xj ) =

τ=

Nd max(Dc )
×
Nc max(Dd )

(20)

Where Nd is the number of discrete features and Nc is the amount of continuous features. max(Dc ) stands for the maximum value of continuous distance
and max(Dd ) demotes the maximum value of discrete features. And the
overall distance will be computed by:
D(xi , xj ) = Dc (xi , xj ) + τ ∗ Dd (xi , xj )
12

(21)

This method is applied in our program, the details will be further illustrated
in section 5.
4. Mathematical Analysis of the Proposed Method
The novelties of our method mostly lie in the designing of the embedded
graph with Gaussian and Delta Kernel (equation 8) and the clustering-based
sampling method (equation 13). The rationality of the new graph-neightbour
method is easy to prove, as the kernel methods used are symmetric and positive semi-definite. However, one might have doubt of our proposed method
and its assumption and rationality. In this section, we will analyze the assumptions and the mathematical properties of the proposed method, and
show that the method could have an analytical upper bound of error rate.
The basic assumption has been stated in section 3.2 with equation 5. Here
to further illustrate the assumption, we could denote the sample number of
underlying patterns and the classes with the following matrix:
 
  
N1
p11 p12 · · · p1n
c1
N2   p21 p22 · · · p2n   c2 
 =
 
(22)
· · · · · · · · ·
· · ·  · · ·
Nd
pd1 pd2 · · · pdn
cn
Where Ni means the sample number of class i = 1, 2, ..., d and cj means
the number of samples belonging to underlying pattern j ∈ 1, 2, ..., n. The
component pij of the d × n matrix P denotes the proportion of underlying
cluster j to be assigned to category (treatment) i. According to equation 10,
data after clustering could also be partitioned into different parts with the
d × n matrix W . Now we assume that the corrupted data (with mislabeled
samples) could be denoted as:
 ∗  ∗
 
N1
p11 p∗12 · · · p∗1n
c1
N2∗   p∗21 p∗22 · · · p∗2n   c2 
 =
 
(23)
 · · ·  · · · · · ·
· · ·  · · ·
cn
Nd∗
p∗d1 p∗d2 · · · p∗dn
And here P ∗ = p∗ij could be regarded as the true distribution proportion
matrix with corrupted (misclassified/mistreated) samples. And our clustering method is performed upon this matrix. To satisfy our assumptions, the
matrix should have the following property:

13

d
X

p∗ij = 1

(24)

i=1

To make sure that all samples from each cluster will be fully assigned to
different classes. And based on our assumption in equation 3, the following
property should hold as the number of the samples asymptotically grows to
infinity:
X
X
(25)
p∗ki ci >
p∗kj cj , ∀k ∈ 1, 2, ..., d
i∈I

j∈I

Where I is the set of the pattern symbols that pki > 0 and J is the set with
originally pkj = 0, which means that cluster j should not have any member
been assigned to class k. For each class k ∈ 1, 2, .., d, the whole number of
data is composed of:
Nk∗

=

X

p∗ki
Pd

i∈I

k=1

p∗ki

+

p∗kj

X

Pd
j∈J

k=1

p∗kj

(26)

In the above equation, the first term of the right hand side comes from the
clusters that should be assigned to the class and the second term comes as
pure mistakes. The error-specified size of sample could then be calculated as:
ε=

X |pki − p∗ |ci X p∗kj cj
ki
+
Pd
Pd
∗
∗
p
ki
k=1
k=1 pkj
j∈J
i∈I

(27)

Plugging in the equation 24 and 25 to equation 27, we could get:
X
X
ε=
|pki − p∗ki |ci +
p∗kj cj
i∈I

<

X

j∈J

|pki −

p∗ki |ci

i∈I

=

X
i∈I1

+

X

p∗ki ci

(28)

i∈I

pki ci +

X

2(p∗ki − pki )ci

i∈I2

Where I1 stands for the set that p∗ki 6 pki , and I2 denotes the set of i that
p∗ki > pki . And by computing the error rate η with ErrorRate = ErrorSamples/(ErrorSamples + CorrectSamples), we could get the following result

14

from equation 27:
P
P
∗
i∈I1 pki ci +
i∈I2 2(pki − pki )ci
P
η<P
∗
∗
i∈I1 (pki ci + pki ci ) +
i∈I2 2pki ci

(29)

And now, we could further analyze equation 29 with two special cases.
1. When p∗ki 6 pki , ∀i ∈ 1, 2, ..., n.
P
Then we will have a number of correct samples as i∈I p∗ki ci . The error
rate would therefore be:
P
i∈I pki ci
(30)
η<P
∗
i∈I (pki ci + pki ci )
And equation 30 would be the upper bound of the error rate under this
situation. As we could see, in this case the function is a decrement
function w.r.t p∗ki , and when p∗ki = pki ∀i, which means at this special
point, the upper bound of error would have a minimal value of η < 0.5.
2. When p∗ki > pki ∀i ∈ 1, 2, ..., n.
P
Then we will have a number of correct samples as i∈I pki ci . The error
rate would therefore be:
P
2(p∗ − pki )ci
P ki ∗
(31)
η < i∈I
i∈I 2pki ci
The equation 31 is the upper bound of the second situation and it is
an increment function. As we could prove the minimal of this upper
bound is also 0.5, and the situation is still p∗ki = pki ∀i. This comes with
the property of a continuous upper bound of the method.
And in the above method we have shown the assumptions of the method and
why it is rational. In practice, if the clustering algorithm is strong enough
to capture the majority of correct-specified samples, then p∗ki should be close
to pki given the assumption in 3.
5. Practical Settings and Experimental Results
In addition to algorithm designing and analyzing, in this work parameter
setting problems are also discussed and experimental results are obtained
based on SEER data. This section will first introduce the SEER data, and
then discuss (hyper)parameter setting problems. And finally, the experimental results with our method will be demonstrated.
15

5.1. SEER Dataset
The SEER dataset from 2004 to 2014 was queried for patients who were
confirmed pancreatic adenocarcinoma defined according to the International
Classification of Disease for Oncology, Third Edition (ICD-O-3) codes for
morphology (8140 and 8500) and topography (C25.0, C25.1, C25.2, C25.3,
C25.7, C25.8, and C25.9). Data collected for each patient includes patient
characteristics (age at diagnosis, sex, race, and marital status), tumor characteristics (tumor location, tumor grade, and AJCC stage), and treatment
characteristics (type of surgery, radiotherapy and status of surgery adjuvant radiotherapy). These three treatments are considered as the three
classes/categories under the term of Machine Learning in out work.
There are some missing data in the dataset, and discussing them will be
out of our research scope. Thus, in our program we only consider the
features(covariates) without data missing, and the rationality of the selection of features could be supported from the perspective of bioscience and
medical science. According to the true meaning of each feature, there only
four features which could be regarded as continuous: EOD10 PN, CS SIZE,
CS EXT, and CS NODE. Some of the features could imply the treatment
(say, NO SURG flag) or should not be treated as feature (say, CASENUM,
the index of patient), so we also remove these features.
The original data is believed to only have mild mis-labeled samples (mistreated patients). To demonstrate the performance of our method under data
with different extent of corruption, we also generate some further-corrupted
data for training purpose. We produce data with 10%, 20% and 40% of artificial label corruptions based on the original data by corrupting the labels
to other two classes with equal opportunities. Notice that the artificially
corrupted data is only utilized for training, and when it comes to evaluation,
we will use the original data to test the robustness of our proposed approach.
5.2. Parameter Settings
Following the method proposed in section 3, there are some hyper- parameters that should be assigned based on our assumption and practical
Requirement. Here we consider three aspects of parameter setting, namely
the number of clusters, Spectral Clustering and Xgboost parameter settings,
and considerations of tSNE.
The first hyper-parameter to determine is the number of clusters we would
like to use. This reflects the our basic assumption of the number of underlying patterns of the data. In our program we choose n = d where n is the
16

number of clusters and d as the number of categories (classes/treatments).
This setting is for a better interpretation property: if there are no overlaps
between each classes (mutually exclusive, say partial gastrectomy and total
gastrectomy), then ideally the clustering result should have a distinct distribution between different classes and the correctly labeled samples will be
straightforward to discover; otherwise, if there are overlaps between different
classes/treatments(say in our SEER data situation, there are operation treatment, radian treatment and operation+radian treatment), then this setting
will be convenient for us to analyze the demanded components of each class.
In the subsection 5.3 we will demonstrate the interpretation of the SEER
data we used in the experiment based on this scheme of cluster number setting.
The second topic of parameter setting is our selections for the Spectral Clustering and Xgboost API. In our work, we use the Spectral Cluster method
provided by SK-learn package with Python API. In the package, we specify
to use Arpack to compute the eigenvalue decomposition. This will lead to
a slightly higher time consumption, but the result will be significantly more
stable. The final clustering results are obtained by K-mean over the embedded space, and to obtain a close-to-optimal result, the method will compute
K-means in the embedded space with 10 different initializations to get the
best result. And in our program, the σ value used in equation 6 could either
be specified or be computed via equation 7. The parameters of Xgboost has
been discussed in section 3.3. In practice, we set tree maximum depth as
7, and η value as 0.5 to control the learning with a relatively conservative
scheme. We specified the max training iteration as 10 rounds for original
data and 5 rounds for artificially corrupted data, as this setting could provide as an optimal validation error.
The third point to be considered in practice is the tSNE visualization technique. In this project, we use the tSNE package provided by SK-learn. The
authors of the original paper proposed tSNE suggested that hyper-parameters
of tSNE is robust and won’t make too large an effect on the result. However,
in our practical experiment, things are much more tricky: the parameter of
’perplexity’, which determines the balance between local manifold and global
manifold, could significantly affect the clustering result and the optimal setting of this parameter will vary with the number of samples for visualization.
Briefly speaking, the quantity of perplexity should be in the same magnitude
of number of samples. For example, if we have 300 samples, then a perplexity
of 5-50 will be good for visualization; however, if we have 30000 sample, then
17

it seems the value of perplexity score should be adjusted to a larger value
accordingly.
In our program, we visualized the project of data in different clusters and
treatments(classes) with tSNE respectively. For the cluster illustration, we
randomly choose 50 samples from each cluster (150 sample in total) and set
perplexity=90, learning rate = 7 and use our precomputed distance with
equation 21; And for the classes illustration problem, we choose 100 samples from each treatment group (before and after processed by our sampling
method) and set perplexity=180 and learning rate = 7.
5.3. Experimental Result
There are in total 8683 valid records, with 3531 being treated with method
1, 3200 with method 2 and 1952 with method 3. And the result of clustering divide the data into a proportion of 2905 in cluster 1, 3817 in cluster 2
and 1961 in cluster 3. Specifically, the treatment-cluster data segment result
could be shown in table 1:
And the sampling probability computed by equation 13 could be illustrated
Cluster 1
New Treatment 1 210
New Treatment 2 2615
New Treatment 3 80

Cluster 2
2159
250
1408

Cluster 3
1162
335
464

Table 1: Treatment-Clustering Result

in table 2.
As we mentioned in the previous sections, the setting of d = n could helps us
Treatment 1
Treatment 2
Treatment 3

Cluster 1
0.0587
0.7921
0.0435

Cluster 2
0.4597
0.0576
0.5827

Cluster 3
0.4816
0.1503
0.3738

Table 2: Probability to sample for new data

explain the result. In table 1 and 2, it could be derived that: 1. patient characterized by cluster 1 should be predominately treated with treatment
2 (radian therapy); 2. patients characterized by cluster 3 and cluster 1
prefer to be treated by a treatment involved with operations, but those
18

who belongs to cluster 3 are more likely to reject radian therapy.
The illustration of the clustering result in 2-dimension could be shown as
figure 1. The result is produced with the settings stated in section 5.2 and
we could find out that the three clusters have clear boundaries against each
other, implying that the result should be rational. The red points, which
represent data from cluster 1, have a more distinct boundary comparing to
cluster 2 and 3, and this could be a supporting evidence for patients characterized by cluster 1 to be treated by a separate method.
The 2-d projection of training sets before and after process could be il-

Figure 1: Clustering Result in 2d. In the figure, red samples are from cluster
1, blue samples are from cluster 2 and green samples are from cluster
3. The group of red samples has a clear boundary with the rest two clusters,
comparing the relatively vague boundary between cluster 2 and 3

lustrated in figure 2. The figure could prove the rationality of our process
and the overlap between treatment 1 and 3. Meanwhile, we could find out
the data manifolds are more compact and have multiple ’clusters’ inside the
class, indicating that the dataset after process could be better categorized
comparing to those before process. This is more significant in figure 3 as
we specifically illustrate the manifold of treatment 1 and treatment 3.
To evaluate the quality of the computed propensity score, we need to compute the confounding metric after applying propensity score to the original
data. Here we consider propensity score weighting and employ Standardized Bias (SB) for several key features under ATE settings according to
[17]. The metric could denote the level of confound between two groups
(treatment and control groups) for each feature (covariate) and if the confounding level is lower (smaller in the metric) then it indicate that the data
after propensity score weighted is better randomized. The evaluation method
19

Figure 2: The comparison of data projection in 2-d plane of training set before
and after the proposed process. left: data manifold before process; right: data
manifold before process. In the figure, red samples are from treatment 1, blue
samples are from treatment 3 and green samples are from treatment 3.

Figure 3: The comparison of data projection inside treatment 1 and 2 in 2-d
plane of training set before and after the proposed process. In the figure, red
samples are from treatment 1 and green samples are from treatment 3.
top: data manifold before process; bottom: data manifold before process. From
the figure, we could find out that data manifold after our process tends to have
some difference between treatment 1 and 3, although they are still very similar.

20

(a) Comparison of data (b) Comparison of data (c) Comparison of data
with 20% corrupted
with 10% corrupted
with 40% corrupted
Figure 4: Data manifold before and after sampling procedure for data that added
artificial label corruptions. red samples are from treatment 1, blue samples
are from treatment 3 and green samples are from treatment 3. For each
sub-figure, left: before precess, right: after process

Age
Processed Data + 0.165
Xgboost
Original Data + 0.138
Xgboost
Raw Data with- 0.089
out
Propensity
Weighting

Race

CS
NODE
0.0025 0.0085

CS
SIZE
0.2054

D
AJCC S
0.0695

Survive
Month
0.3118

0.0241 0.1096

0.0519

0.1377

0.1315

0.0446 0.2346

0.3857

0.7434

0.5806

Table 3: Standardized Bias between treatment group 1 and 2


could only compare two groups (classes) so that here we have 32 = 3 combinations. The comparison of the Standardized Bias could be found in table
3, 4 and 5. Here, the Processed Data was trained with xgboost for 10 iterations and the Original Data is a random subset of 3000 samples from the
whole dataset. Final training and testing mlogloss on Processed Data are
0.243960 and 0.552677; Final training and testing mlogloss on Original Data
are 0.222134 and 0.500497.
To demonstrate the effect of noise-correction for further corrupted data, we
also compared the Standardized Bias between treatment 1 and treatment 2
for the data with 10%, 20% and 40% artificial corrupted label in table 6,
7 and 8 respectively. Notice here when computing the Standardized Bias
with the whole dataset, we use the authentic data. This will be the common

21

Age
Processed Data + 0.052
Xgboost
Original Data + 0.134
Xgboost
Raw Data with- 0.137
out
Propensity
Weighting

Race

CS
NODE
0.0085 0.1597

CS
SIZE
0.2363

D
AJCC S
0.1560

Survive
Month
0.0218

0.0051 0.3821

0.4261

0.3287

0.0714

0.0605 0.0364

0.0875

0.0307

0.1113

Table 4: Standardized Bias between treatment group 1 and 3

Age
Processed Data + 0.238
Xgboost
Original Data + 0.297
Xgboost
Raw Data with- 0.246
out
Propensity
Weighting

Race

CS
NODE
0.0134 0.0981

CS
SIZE
0.3440

D
AJCC S
0.2042

Survive
Month
0.3237

0.0357 0.1130

0.3157

0.4209

0.2432

0.0195 0.1997

0.4235

0.8343

0.8037

Table 5: Standardized Bias between treatment group 2 and 3

setting in real-life situations. Also notice that for the corrupted experiments
the xgboost only computed for 5 iteration for the overfitting consideration.
From the table, we could find that 1. propensity score computed by processed
data with xgboost could usually outperform the scores computed by original
sampled data with the same method; 2. the advantages of processed data will
increase as the proportion of corrupted data increase 3. xgboost generally
has a excellent and robust performance in computing propensity scores, and
the propensity score computed could still improve the randomization even
with 40% corrupted labels.
6. Conclusion
In this paper, a novel approach to deal with the corrupted-label problem
in computing propensity score for medical data is proposed and its property
22

Age
Processed Data +
Xgboost
Original Data +
Xgboost

Race

0.053

CS
NODE
0.0269 0.0932

CS
SIZE
0.3351

D
AJCC S
0.4977

Survive
Month
0.3784

0.045

0.0666 0.0811

0.1496

0.4988

0.4201

Table 6: Standardized Bias between treatment group 1 and 2 with 10% labels
corrupted

Age
Processed Data +
Xgboost
Original Data +
Xgboost

Race

0.059

CS
NODE
0.0553 0.0149

CS
SIZE
0.2696

D
AJCC S
0.4712

Survive
Month
0.4745

0.037

0.0474 0.1238

0.2757

0.5243

0.5577

Table 7: Standardized Bias between treatment group 1 and 2 with 20% labels
corrupted

is analyzed. The Xgboost approach is utilized to compute the propensity
score and the performance of the integrated approach is examined and compared. TSNE method is employed for comprehensive purpose and the data
manifold of different clusters and treatments are demonstrated.
The paper made following major contributions to the field: Firstly, the paper finds out the mislabeling-data problem in medical science and researched
on this topic. Given the fact that the problem has been under-investigated
before but the phenomenon is very common, the research would significantly
fill the blank of this aspect in bioscience research; Secondly, the paper proposed an effective method for getting a relatively authentic subset of training
data, and proved mathematically about its rationality and effectiveness. The
rigorous proof makes it possible for this method to be utilized into general
Machine Learning problems and solve a broader scope of tasks. Thirdly, the
paper employed Xgboost and tSNE methods to this area, which are relatively
novel attempts. The results imply a huge research potential in this area.
In the future, the author would like to concentrate more on mislabeling data
in bioscience area. We intend to try penalization-based label correction methods proposed in the field of Machine Learning, and design new algorithms
23

Age
Processed Data +
Xgboost
Original Data +
Xgboost

Race

CS
NODE
0.0250 0.0197 0.0968

CS
SIZE
0.3667

D
AJCC S
0.5715

Survive
Month
0.4491

0.0968 0.0244 0.1230

0.3194

0.5793

0.5826

Table 8: Standardized Bias between treatment group 1 and 2 with 40% labels
corrupted

specifically in dealing with medical data. Meanwhile, to analyze the confounding covariates, we plan to research on graphical model-based methods
in the long run.

24

References
[1] P. R. Rosenbaum, D. B. Rubin, The central role of the propensity score
in observational studies for causal effects, Biometrika 70 (1983) 41–55.
[2] W. Meng, S. Wang, J. Xia, J. Lyu, F. Shi, A simulation study for the
treatment effect using propensity score matching, Unpublished (2015).
[3] T. Chen, C. Guestrin, Xgboost: A scalable tree boosting system, CoRR
abs/1603.02754 (2016).
[4] Propensity score estimation: neural networks, support vector machines,
decision trees (cart), and meta-classifiers as alternatives to logistic regression, Journal of Clinical Epidemiology 63 (2010) 826 – 833.
[5] J. Kruppa, Y. Liu, G. Biau, M. Kohler, I. R. Knig, J. D. Malley,
A. Ziegler, Probability estimation with machine learning methods for
dichotomous and multicategory outcome: Theory, Biometrical Journal
56 (2014).
[6] A. Linden, P. R. Yarnold, Using classification tree analysis to generate
propensity score weights, Journal of Evaluation in Clinical Practice 23
(2017) 703–712. JECP-2017-0075.
[7] B. K. Lee, J. Lessler, E. A. Stuart, Improving propensity score weighting
using machine learning, Statistics in medicine 29 (2010).
[8] van der Laan, E. M., Polley, A. Hubbard, Super learner, Statistical
Applications in Genetics and Molecular Biology 6 (2007).
[9] R. Pirracchio, M. L. Petersen, M. van der Laan, Improving propensity score estimators’ robustness to model misspecification using super
learner, American Journal of Epidemiology 181 (2015) 108–119.
[10] N. Kreif, R. Grieve, I. Daz, D. Harrison, Evaluation of the effect of a
continuous treatment: A machine learning approach with an application
to treatment for traumatic brain injury, Health Economics 24 (2015)
1213–1228.
[11] C. Ju, M. Combs, S. D. Lendle, J. M. Franklin, R. Wyss, S. Schneeweiss,
M. J. van der Laan, Propensity score prediction for electronic healthcare
25

databases using Super Learner and High-dimensional Propensity Score
Methods, ArXiv e-prints (2017).
[12] C. Leyrat, S. R. Seaman, I. R. White, I. Douglas, L. Smeeth, J. Kim,
M. Resche-Rigon, J. R. Carpenter, E. J. Williamson, Propensity score
analysis with partially observed confounders: how should multiple imputation be used?, ArXiv e-prints (2016).
[13] R. Mitra, J. P. Reiter, Estimating propensity scores with missing covariate data using general location mixture models, Statistics in Medicine
(2011) 627–41.
[14] N. Natarajan, I. S. Dhillon, P. K. Ravikumar, A. Tewari, Learning with
noisy labels, in: C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani,
K. Q. Weinberger (Eds.), Advances in Neural Information Processing
Systems 26, Curran Associates, Inc., 2013, pp. 1196–1204.
[15] Y. Li, J. Yang, Y. Song, L. Cao, J. Luo, J. Li, Learning from noisy
labels with distillation, CoRR abs/1703.02391 (2017).
[16] A. de Brébisson, P. Vincent, An exploration of softmax alternatives
belonging to the spherical loss family, CoRR abs/1511.05042 (2015).
[17] D. F. McCaffrey, B. Griffin, D. Almirall, M. Slaughter, R. Ramchand,
L. Burgette, A tutorial on propensity score estimation for multiple
treatments using generalized boosted models, Statistics in medicine 32
(2013) 3388–3414.

26

