Dealing with the Data Imbalance Problem in Pulsar Candidate Sifting Based on Feature Selection

A pulsar candidate is originally a piece of signal from the receiver of a radio telescope during the observation time. Most commonly, it is processed by PulsaR Exploration and Search TOolkit (PRESTO; Ransom 2001), a typical software for pulsar search and analysis. Then the candidate is represented using a series of physical values and a series of diagnostic plots as Figure 1 shows. On the left, the plots from top to bottom are: a subband plot, a sub-integration plot and a folded profile of the signal. A subband plot displays the pulse in different bands of observed frequencies; a sub-integration plot shows the pulse in the time domain; a folded profile is the folded signal of its subbands by frequency or the folded signal of its sub-integrations by period. On the right are two plots. One is a grid searching plot for DM and period. The other is a DM-searching curve. It is known that DM measures the number of electrons which a pulsar's signal travels through from the source to the Earth. However, a real DM is unknown and should be obtained by trials. Then a grid searching plot on the right top of Figure 1 records the change of signal-to-noise ratio (S/N) as the trial DMs and the trial periods vary. A DM searching curve on the right middle describes the relationship between trial DM and its corresponding S/N, and the peak of the curve implies the most likely value of DM.

Our work is conducted on the HTRU survey. The HTRU data set is observed with two observatories, the Parkes radio telescope in Australia and the Effelsberg 7-beam system in Germany. HTRU is an ambitious (6000 hr) project to search for pulsars and fast transients in the entire sky, which is split into three areas: Low-latitude (covers ±35), Mid-latitude (Medlat; covers ±15°) and High-latitude (covers the remaining sky <10°). The pipeline searched for pulsar signals with DMs from 0 to 400 pc·cm⁻³ (DM is often quoted in the units of pc·cm⁻³, which makes it easy to estimate the distance between a given pulsar and the Earth), and also performed an accelerated search between −50 to +50 m s⁻².

The HTRU data set was publicly released by Morello et al. (2014) and is available online ³ . It consists of 1196 real pulsar candidates and 89,995 non-pulsars, which are highly class-imbalanced, as only a tiny fraction of the candidates are from real pulsar signals (Table 1).

PCS can be described as a binary supervised classification issue in ML. Supervised learning (Mitchell et al. 1997) is an ML task of learning a function that maps instances to their labels. Particularly, a classifier of PCS aims to learn function mapping features of the pulsar candidates to their categories—-pulsar or non-pulsar. To evaluate the effectiveness of selected features, the performance of classifiers should be estimated.

In our work, feature selection algorithms are evaluated by seven classifiers. Among them, DT (Quinlan 2014), LR (Hosmer et al. 2013) and SVM (Suykens & Vandewalle 1999) are normal classifiers, while Adaptive boosting (Adaboost; Freund & Schapire 1997), GBC (Möller et al. 2016), XGboost (Chen & Guestrin 2016) and Random Forest (RF; Liaw et al. 2002) are ensemble learning classifiers. The principles of these classifiers are different and representative. For example, a typical DT classifier is based on the information gain ratio while SVM tries to find the best hyperplane which represents the largest separation between the two classes. Ensemble methods (Dietterich et al. 2002) use multiple weak classifiers such as DT to obtain a strong classifier. Interested readers can refer to their references for details (Mitchell et al. 1997; Mohri et al. 2018).

To evaluate the performance of a classifier on class-imbalanced data, typically on pulsar candidate data, four most relevant metrics are given. They are the Recall rate, Precision rate, F₁ score and FPR. They can be expressed by True Positives (TP), True Negatives (TN), False Positives (FP) and False Negatives (FN).

In binary classification, recall, defined by TP/(TP + FN), where TP denotes the amount of true pulsars predicted and TP + FN the total amount of true pulsars. It measures how many pulsars could be correctly predicted pulsars from all the real pulsars. Precision, defined by TP/(TP + FP), measures how many true pulsars would be predicted correctly out of the candidate prediction tasks. However, recall and precision compose a relationship that opposes to each other, as it is probable to increase one at the cost of reducing the other. Therefore, F₁ score, defined as the harmonic mean of recall and precision, i.e., (2 · Precision · Recall)/(Precision + Recall), is a trade-off between them. As for FPR, it measures the ratio of mislabelled non-pulsars out of all the non-pulsar candidates by FP/(TN + FP). It can be inferred by recall and precision. Therefore, we just need to focus on the recall, precision and F₁ score of a classifier.

To extract some robust and useful features, guidelines of feature design were proposed by PCS researchers. Morello et al. (2014) gave several suggestions, such as "ensuring complete robustness to noisy data" in order to "exploit properly in the low-S/N regime," and "limiting the number of features" to deny "the curse of dimensions." Lyon et al. (2016) suggested that features should be designed to maximize the separation between positive and negative candidates, reducing the impact of class imbalance.

Based on suggestions from Morello et al. (2014) and Lyon et al. (2016), guidelines for candidate features in a feature pool were summarized as follows. Candidate features:

• 1.  
  Should be distinguishable enough between pulsars and non-pulsars. A distinguishable feature will greatly improve the performance of the classifier.
• 2.  
  Should be diversified. Considering the diversity of the feature source, both empirical features and statistical features should be included in the feature pool.
• 3.  
  Should fully cover the space. Features should be extracted from the mainly diagnostic images, especially the sub-integration plot, the subband plot, the folded profile and the DM searching curve.
• 4.  
  Can be easily extracted and calculated.
• 5.  
  Should be controlled in a moderate total number. If the total number is small, some relevant features may be missed; if it is large, it will enlarge the computing cost of the feature selection algorithms.

Following the guidelines in Section 3.1, 22 features have been collected (Table 2), which are candidate features from Morello et al. (2014), Lyon et al. (2016) and Tan et al. (2017). Among them, six features were defined by Morello et al. (2014) to build the SPINN model, eight statistical features were proposed by Lyon et al. (2016) as inputs to GHVFDT and eight additional features were introduced by Tan et al. (2017) to develop an ensemble classifier comprised of five different DTs. Details of these features are described in Table 2.

Table 2. Notations and Definitions of 22 Candidate Features in Our Work

ID	Feature	Description
M1	log(S/N)	Log of the signal-to-noise of the pulse profile
M2	Deq	Intrinsic equivalent duty cycle of the pulse profile
M3	$\mathrm{Log}(P/\mathrm{DM})$	Log of the ratio between period and DM
M4	VDM	Validity of optimized DM
M5	χ(S/N)	Persistence of the signal in the time domain
M6	Drms	rms between folded profile and sub-integration
L1	Pfμ	Mean of the folded profile
L2	Pfσ	Standard deviation of the folded profile
L3	Pfk	Kurtosis of the folded profile
L4	Pfs	Skewness of the folded profile
L5	DMμ	Mean of the DM curve
L6	DMσ	Standard deviation of DM curve
L7	DMk	Kurtosis of DM curve
L8	DMs	Skewness of DM curve
T1	SubbandCorrμ	Mean of the SubbandCorr1
T2	SubbandCorrσ	Standard deviation of SubbandCorr
T3	SubbandCorrk	Kurtosis of SubbandCorr
T4	SubbandCorrs	Skewness of SubbandCorr
T5	SubintCorrμ	Mean of the SubintCorr2
T6	SubintCorrσ	Standard deviation of SubintCorr
T7	SubintCorrk	Kurtosis of SubintCorr
T8	SubintCorrs	Skewness of SubintCorr

Note. Features with ID M1–M6 were defined by Morello et al. (2014); Features with ID L1–L8 were created by Lyon et al. (2016); Features with ID T1-T8 were defined by Tan et al. (2017). 1. SubbandCorr: A vector of correlation coefficients between subband and the folded profile. 2. SubintCorr: A vector of correlation coefficients between sub-integration and the folded profile.

Download table as:  ASCIITypeset image

To demonstrate the discriminating capabilities of these features, one statistical approach is to show the distributions of pulsars and non-pulsars from each feature by box plots, which can graphically demonstrate the locality, spread and skewness groups of the features. Figure 2 displays box plots of our candidate features on HTRU. There are two box plots per feature. The red boxes describe the feature distribution for known pulsars, while the white ones are for non-pulsars mainly consisting of RFI. Note that the data of each feature were scaled by z-score, with mean zero and standard deviation one. The resulting z-score measures the number of standard deviations that a given data point is from the mean. Generally, the less the overlap of the red box and white box in a feature, the better the separability of the feature. However, the usefulness of features according to their box plots is only on a visual level. Measurable investigation of these features will be given in the next section.

Zoom In Zoom Out Reset image size

A feature selection algorithm is a search technique for a feature subset from a feature pool. Irrelevant or redundant features not only decrease the calculation efficiency of ML models but also negatively affect their performances. Better selected features can be helpful when facing a data-imbalanced problem as described by Wasikowski & Chen (2010), and some practical algorithms have been proposed for a two-class imbalanced data problem (Yin et al. 2013; Maldonado et al. 2014).

There are main three categories of feature selection algorithms: filters, wrappers and embedded methods. Filter methods use a proxy measure to score a feature subset. Commonly, they include the mutual information (Shannon 1948), the point-biserial correlation coefficient (Gupta 1960) and Relief (Urbanowicz et al. 2018) score. Wrapper methods score feature subsets using a predictive model. Common methods include grid search, greedy (Black 2012) and recursive feature elimination (Kira & Rendell 1992). An embedded feature selection method is an ML algorithm that returns a model using a limited number of features.

Our proposed algorithm of feature selection will combine the idea of a filter called Relief with a wrapper method named Greedy. On one hand, PCS, as a mission of binary supervised classification, is highly dependent on the constructive features between pulsars and non-pulsars. Thus, methods of filtering are preferably considered as they measure the relation between features and labels. In fact, measures by filters are able to capture the usefulness of the feature subset only based on the data, which are independent of any classifier. A Relief algorithm is one of the best measures of filters when compared with other filter methods. It weights features and avoids the problem of high computational cost in a combinatorial search. Thus, our proposed approach of feature selection for PCS is a Relief-based algorithm. On the other hand, the selected features according to their Relief scores may be redundant. They involve more than two features with high Relief scores but they are strongly correlated, since one relevant feature may be redundant in the presence of another relevant feature (Guyon & Elisseeff 2003). To remove these redundant features, it follows the evolution of Relief scores to implement a Greedy technique.

However, Greedy or Relief for feature selection has some shortcomings. Although Relief score is able to filter out some irrelevant features, it could not detect the redundant ones, which implies that if two features share the same information in terms of correlation measure, both of them are very likely judged as relevant or irrelevant. As for Greedy, it can be utilized to reduce the number of features. Its computational cost increases in a quadratic way as the number of features increases, which makes the computation unaffordable.

Based on the considerations above, we combine Relief with the Greedy algorithm to propose a KFRG algorithm of feature selection for PCS. In the first stage, it filters out some irrelevant features according to their Relief scores, while in the second stage, it removes the redundant features and selects the most relevant features by a forward Greedy search strategy. Experimental investigations on HTRU showed that it improved the performance of most classifiers and achieved high both recall rate and precision (Section 5).

Relief (Kira & Rendell 1992) is a filter algorithm of feature selection which is notably sensitive to feature interactions. It calculates a feature score for each feature which can then be applied to rank and select top scoring features for feature selection. Alternatively, these scores may be applied as feature weights to guide downstream modeling. Relief is able to detect conditional dependencies between features and their labels (pulsars and non-pulsars) and provide a unified view on the feature estimation in regression and classification. It is described as Formula (1). The greater the Relief score of a feature is, the more distinguishable the feature is, and corresponding features are more likely to be selected.

Let $D=\{({X}_{i},{y}_{i})| {X}_{i}=({x}_{i}^{1},{x}_{i}^{2},\cdots ,{x}_{i}^{d}),i\in I\}$ denote a data set, where y_i is the label of X_i and x_i ^j the jth component of X_i . Denote ${x}_{i,{nh}}^{j}$ as the jth feature of its nearest instance whose label is the same as that of X_i (a "nearest hit"), while ${x}_{i,{nm}}^{j}$ the jth feature of its nearest instance whose label is different from that of X_i (a "nearest miss"). Then Relief score of the jth feature (denoted as δ^j ) is defined as

$$\begin{eqnarray}&&{\delta }^{j}=\displaystyle \frac{1}{| I| }\displaystyle \sum _{i\in I}(-{diff}{({x}_{i}^{j},{x}_{i,{nh}}^{j})}^{2}+{diff}{({x}_{i}^{j},{x}_{i,{nm}}^{j})}^{2}),\end{eqnarray} \tag{ 1 }$$

where diff represents the difference of two components. ${diff}({x}_{a}^{j},{x}_{b}^{j})=\left\{\begin{array}{ll}0, & {y}_{a}={y}_{b}\\ 1, & {y}_{a}\ne {y}_{b}\end{array}\right.$ if the jth feature is discrete, while ${diff}({x}_{a}^{j},{x}_{b}^{j})=| {x}_{a}^{j}-{x}_{b}^{j}|$ if the jth feature is continuous.

Algorithm 1. Relief Algorithm

Input:

D; #Training data set

S; #Features pool

δ0; #A preset threshold

Output:

Δ. Set of selected features

1: ${\rm{\Delta }}=[\ ]$ # Initialization of the chosen features

2: for $j=1:| S|$ do

3: Compute ${\delta }^{j}$ by Equation (1)

4: if ${\delta }^{j}\gt {\delta }_{0}$ then

5: ${\rm{\Delta }}={\rm{\Delta }}\cup \{j\}$

6: return Δ

Download table as:  ASCIITypeset image

A greedy algorithm is an algorithm that follows the problem-solving heuristic of making the locally optimal choice at each stage (Black 2012). The greedy strategy does not usually produce an optimal solution. Nonetheless, a greedy heuristic method may yield a locally optimal solution that approximates the global optimal solution in a reasonable amount of time.

In our algorithm, the objective function (target) is to maximize the F₁ score of a classifier, and thus our greedy implementation is designed to choose the best feature step by step from the rest of the feature pool. Here, "the best feature candidate" is the feature which contributes best to increase of the F₁ score of the classifier. Accordingly, the stopping criterion of our greedy implementation in the iteration procedures is that the F₁ score does not increase anymore, or the maximum number of selected features is more than a preset threshold Maxlen.

Maxlen is a hyperparameter to control the maximum number of selected features. On one hand, Maxlen should be large enough to enable as many candidate features as possible. A small Maxlen may miss some relative features and result in low performance. On the other hand, as Maxlen increases, so does the computational complexity (Goldreich 2008). Computational complexity is an important part of an algorithm design, as it gives useful information about the amount of resources required to run it. In fact, the computational complexity of Greedy can be expressed as $O({{Maxlen}}^{2})\times O({\mathfrak{L}})$ according to Algorithm 2, where the notation O(Maxlen²) means the run time or space requirements grow as the square of Maxlen grows, while $O({\mathfrak{L}})$ represents the computational complexity of ${\mathfrak{L}}$ which relies on both the choice of a classifier and the size of the input (feature). Thus, a large Maxlen implies a large computing cost. To evaluate a fit, Maxlen experiments with values ranging from 2 to 10 were carried out. Experimental investigation shows that as Maxlen increases, the average performance metrics improve rapidly at first, while they stay at a similar level when Maxlen is more than 8. The average recall and precision stay around 97.2% and 98.0% respectively for Maxlen larger than 8 as Table 5 affirms. Considering both the computational complexity and the performance, Maxlen is set to be 8.

Following the description above, our Greedy Algorithm of feature selection can be described as Algorithm 2.

Algorithm 2. Greedy Algorithm

Input:

D; # Training data set

S; # Features pool

${\mathfrak{L}};$ # Classifier

Maxlen; # Maximum number of selected features

Output:

Sch. # Set of selected features

1: ${S}_{\mathrm{ch}}=[\ ];$ # Initialization of the chosen features

2: ${S}_{\mathrm{un}}=S;$ # Initialization of the unchosen features

3: ${F}_{1}=0;$ # Initialization value of F1 score

4: while $| {S}_{\mathrm{ch}}| \lt {MaxLen}$ do

5: for $j\in {S}_{\mathrm{un}}$ do

6: ${{\mathfrak{L}}}_{j}={\mathfrak{L}}({S}_{\mathrm{ch}}\cup \{j\})$ # add j to a classifier

7: ${j}^{* }={\mathrm{argmax}}_{j\in {S}_{\mathrm{un}}}{F}_{1}({{\mathfrak{L}}}_{j})$ # choose the best feature

8: ${F}_{1}^{* }={F}_{1}({\mathfrak{L}}({S}_{\mathrm{ch}}\cup \{{j}^{* }\})$

9: if ${F}_{1}^{* }\gt {F}_{1}$

10: ${F}_{1}={F}_{1}^{* }$

11: ${S}_{\mathrm{un}}={S}_{\mathrm{un}}\setminus \{{j}^{* }\}$

12: ${S}_{\mathrm{ch}}={S}_{\mathrm{ch}}\cup \{{j}^{* }\}$ # add j to a classifier

13: else

14: break

15: return Sch

Download table as:  ASCIITypeset image

Relief is easy to operate, but it is not very satisfying in some class-imbalanced scenarios, since it may underestimate those features with high discriminative ability that are in the minority, and ignores the sparse distributional property of minority class samples (Yuanyu et al. 2019). That is, a feature with high Relief score pays more attention to the non-pulsars which are in the majority class and it is hard to identify the pulsars which are in the minority class and thus more promising pulsars will be missed.

To overcome these flaws, the K-fold Relief (KFR) algorithm was designed. The key improvement of KFR is to balance the data by recycling the minority and sampling from the majority. First, split the training data into minority samples and majority ones. Then, produce K disjoint subsets from the majority samples randomly, and merge each subset with the minority into K new data sets, each of which is relatively balanced with the same minority samples. Here, K is a preset integer which is normally the ratio of the majority classes to minority classes in the training data set. Finally, calculate the mean of Relief scores of each set. KFR is able to promote the importance of the minority classes for the estimation of relevant features.

Combining KFR with the Greedy algorithm, we get KFRG. KFRG is a two-stage algorithm: the first stage aims to remove some irrelevant features from the candidate features according to their Relief score, while the second stage is designed to select the most relevant features in a greedy way. It can be described in Algorithm 3.

Algorithm 3. KFRG Algorithm

Input:

D; #Training data

S; #Feature pool

${\mathfrak{L}};$ #A classifier

Maxlen; #Maximum number of selected features

Output:

S*. #Set of selected features

1: Split the training data D into minority class Mi (pulsars) and majority class Ma (non-pulsars);

2: Divide Ma into K disjoint subsets, ${Ma}={\cup }_{k=1}^{K}{({Ma})}_{k};$

3: for $k=1:K$ do

4: Obtain the Relief scores of S on the data set ${({Ma})}_{k}\cup {Mi}$ by Algorithm 1 (Relief).

5: Calculate the mean of Relief scores of S to choose a feature subset Δ.

6: Implement Algorithm 2 (Greedy) on Δ by setting the maximum number of selected features as Maxlen.

7: return S*

Download table as:  ASCIITypeset image

5.1.1. Selected Features Based on KFR Scores

The KFR algorithm is the first stage of Algorithm 3, which outputs the mean of Relief scores of the K-fold training sets. The Relief scores standing for the weights of features were calculated and are shown by a bar graph in Figure 3 for HTRU, where both the scores and their ranks are given. To select the more relevant features, a preset threshold will keep the features with higher scores. Here, a ratio of 0.618 is preset to ensure that the number of selected features is more than half of the total number of candidate features. That is, about 12 features out of 22 are considered to be very relevant (blue bars) and 10 others (gray bars) are less relevant.

Zoom In Zoom Out Reset image size

5.1.2. Selected Features Based on KFRG

KFRG feature selection algorithm is an improvement of the Relief algorithm by combining KFR with a greedy algorithm. To demonstrate the effectiveness of KFRG, an ablation study of stack mode is given. Step by step, the performance metrics of all the classifiers were calculated from Relief to KFR, and finally, to KFRG.

Table 4 gives the numerical calculation of recall, precision, F₁ score and FPR for each classifier with three different feature select algorithms—-Relief, KFR and KFRG, while Figure 4 plots their averaged performance metrics. It shows that KFR performs better than Relief, and KFRG performs best of all. For one thing, KFR has better recall rates, better precision and lower FPRs for all classifiers than the original Relief techniques. For example, the recall raises from 94.8% to 96.1%. These improvements come from the step of K-fold Relief operation, as KFR is designed for the imbalance problem. For another aspect, KFRG keeps a precision rate as high as KFR, and raises the recall rate by 0.9%. Further, KFRG achieves a best F₁ score of 97.5%. These improvements come from the step of Greedy as it aims to maximize the F₁ score by removing the redundant features.

Zoom In Zoom Out Reset image size

To verify the effectiveness of KFRG, comparative experiments with different feature groups were carried out, where three subsets of the feature pool were considered as the inputs of the classifiers, including features from Morello et al. (2014), features from Lyon et al. (2016) and features from KFRG. The performance metrics of recall (Rec), precision (Pre), F₁ score and FPR are expressed in Table 5.

The experimental results in Table 5 are summarized as follows.

• 1.  
  Both the recall and precision have significantly improved.The average recall of the classifiers is 97.2%, and most of the classifiers achieve recall rates ranging from 96% to 99% after applying the KFRG algorithm, which implies that most of the real pulsar signals were well detected after feature selection. For instance, the recall rate and the precision rate of DT are respectively 89.7% and 95.4% with features M1–M6, while they are up to 96.4% and 98% based on KFRG. Also, the average precision is as large as 98.0% based on KFRG, which has increased by an average of 3.9% compared with M1–M6 and L1–L8.
• 2.  
  The FPR of the classifiers is reduced to 0.05% on average.Most of the classifiers achieve an FPR less than 0.05%. A low FPR of a classifier implies that the selected features are very sensitive in excluding non-pulsars.
• 3.  
  F1 scores based on selected features have increased.The best F₁ is 98.3% in our experiment in the following cases: using five selected features [M3, M5, L3, L4, T1], GBDT achieved a recall of 97.8% and a precision of 98.9%; Using another five selected features [M5, M6, L4, T1, T5], XGBoost also achieved a high F₁ score of 98.3%, which is as good as the GBDT classifier. A better F₁ score implies that both recall and precision increase since F₁ score is the harmonic mean of recall and precision. In other words, more potential pulsar signals are correctly recognized and fewer non-pulsar signals are misjudged in these cases.

As feature selection of KFRG alleviates the imbalance problem on ML, the performances based on KFGR were compared with some other widely used data-balancing techniques of oversampling, including random oversampling, SMOTE (Chawla et al. 2002), Borderline SMOTE (Han et al. 2005) and ADASYN (He et al. 2008). Furthermore, we even implement the KFRG-SMOTE method which is a combination of the proposed KFRG and SMOTE techniques. We evaluated the metrics of each algorithm on the different classifiers, and then took the mean of the performance of each classifier on each evaluation statistic to compare between feature selection metrics and oversampling metrics (Table 6).

It shows that KFRG offers better performance than oversampling techniques such as SMOTE, Borderline SMOTE and ADASYN, and random oversampling performs worst among them. KFRG achieves a similar recall rate but an improved precision and a lower FPR, which imply that the classifiers on KFRG are very strict in the criteria for classifying the candidates as pulsars and only a few of the non-pulsars were misjudged. Moveover, the performance of KFRG-SMOTE is as good as KFGR, which shares a high F₁ score of 97.5% and whose FPRs both range at a low level between 0.03% and 0.04%.

Compared with oversampling techniques, KFRG has its own advantages. One of the advantages is that KFRG has good generalization ability and is able to avoid overfitting problems. In fact, random oversampling tends to suffer from overfitting problems as the minority of samples in its training set were duplicated at random. As a result, the trained model becomes too specific to the training data and may not generalize well to new data. Other oversampling techniques are all based on SMOTE, which eliminates the harms of a skewed distribution by creating new minority class samples. They generate a synthetic sample x_new by using the linear interpolation of x and y with the expression of x_new = x + (y − x) × α, where α is a random number in the range [0, 1] and y is a k-nearest neighbor of x in the minority set. However, in many cases the kNN-based approach may generate wrong minority class samples as the above equation says that x_new will lie in the line segment between x and y. In addition, it would be difficult to find an appropriate value of k for kNN a priori. The parameter k varies with the distribution of samples between minority and majority. Different from data-level methods of oversampling, feature selection of KFGR neither duplicates nor generates any additional data. It keeps the distribution and class IR of the data.

The other advantage is that the result of KFRG is interpretable. The selected features are the most distinguishing ones between pulsars and non-pulsars. In addition, KFRG performs well based on the data characteristics regardless of the classifier used. Although the sets of selected features by KFRG may be changed with different classifiers, most of the selected features are the same, as will be explained in the next section.