AUTOMATED CLASSIFICATION OF PERIODIC VARIABLE STARS DETECTED BY THE WIDE-FIELD INFRARED SURVEY EXPLORER

In order to calibrate and validate a classification method for WISE light curves, we assembled a "truth" list of variable stars that were previously classified with measured periods from a number of optical variability surveys. This list includes all eclipsing binaries and RR Lyr stars in the General Catalogue of Variable Stars (GCVS; Samus 2009), the MACHO Variable Star Database (Alcock 2003), and the All-Sky Automated Survey (ASAS; Pojmański 2006). This list of known variables was then positionally cross-matched to the WISE All-Sky Source Catalog (Cutri et al. 2012). A match-radius tolerance of 25 was used and the brightest WISE source was selected if multiple matches were present. Furthermore, only sources with a WISE catalog variability flag (Cutri et al. 2012, Section IV.4.c.iii) of $var\_flg\ge 6$ in the W1 band were retained. This criterion ensured that the WISE source had a relatively good time-averaged photometric signal-to-noise ratio (with typically S/N ≳ 6 in the W1 single exposures) and a high likelihood of also being variable in the mid-IR. After cross-matching, 1320 objects were rejected because they were labeled as either ambiguous (assigned to two or more classes), uncertain, or had different classifications between the optical catalogs. A further 2234 objects had duplicate entries among the optical catalogs (assigned to the same class) and one member of each pair was retained. In the end, we are left with a training sample of 8273 known variables with WISE photometry. Of these 8273 objects, 1736 are RR Lyr, 3598 are Algol-type eclipsing binaries, and 2939 are W UMa-type eclipsing binaries according to classifications reported in previous optical variability surveys.

We assume here that the bulk of classifications reported in the ASAS, GCVS, and MACHO surveys (primarily those objects labeled as Algol, RR Lyr, or W UMa) are correct, or accurate enough as to not significantly affect the overall core-class definitions in the WISE parameter (feature) space. The removal of objects with ambiguous and uncertain classifications, or with discrepant classifications between the catalogs, is expected to have eliminated the bulk of this uncertainty. A comparison of our classification performance to other studies that also used classifiers trained on similar optical catalogs (Section 5.6) shows that possible errors in the classification labels of the training sample are not a major source of uncertainty. However, these errors cannot be ignored.

There is no guarantee that the training sample described here represents an unbiased sample of all the variable types that WISE can recover or discover down to fainter flux levels and lower S/Ns over the entire sky. This training sample will be used to construct an initial classifier to first assess general classification performance (Section 5.6). In the future, this classifier will be used to initially classify WISE flux variables into the three broad variability classes defined above. Inevitably, many objects will remain unclassified in this first pass. To mitigate training sample biases and allow more of the WISE feature space to be mapped and characterized, we will use an active-learning framework to iteratively refine the training sample as objects are accumulated and (re)classified. This will allow more classes to be defined as well as sharpen the boundaries of existing ones, hence alleviating possible errors in the input classification labels (see above). Details are discussed under our general classification plan in Section 6.

A requirement of any classification system is a link between the features used as input and the classes defined from them. We review here the seven mid-IR light-curve features that we found work best at discriminating between the three broad classes that WISE is most sensitive to (defined in Section 2). These were motivated by previous variable star classification studies (e.g., Kinemuchi et al. 2006; Deb & Singh 2009; Richards et al. 2011)

We first extracted the time-resolved mid-IR photometry and constructed light curves for all sources in our training sample from the WISE All-Sky, Three-Band, and Post-Cryo Single-Exposure Source Databases.³ Light curves were constructed using only the W1 band. This is because W1 generally has better sensitivity than W2 for the variable stars under consideration. For each source, a period was estimated using the generalized Lomb–Scargle (GLS) periodogram (Zechmeister & Kürster 2009; Scargle 1982, see Section 4 for this choice). The light curves were phased to the highest peak in the periodogram that fell in the period range 0.126 day (∼3 hr) to 10 days. The lower value corresponds to the characteristic WISE single-exposure sampling and the upper limit is based on the known periods of Algol-type binaries that could be recovered by WISE given the typical time span of observations, with a generous margin. These recovered periods (P_r) constitute our first feature to assist with classification. Section 4 compares these estimates to the available periods derived from optical light curves.

The second feature derived from the WISE light curves in our training sample is the Stetson-L variability index (Stetson 1996; Kim et al. 2011). This index quantifies the degree of synchronous variability between two bands. Because the band W1 and W2 measurements generally have the highest S/N for objects in the classes of interest, only these bands are used in the calculation of this index. The Stetson-L index is the scaled product of the Stetson-J and -K indices:

$$\begin{equation} L = \frac{JK}{0.798}. \end{equation} \tag{ 1 }$$

Stetson-J is a measure of the correlation between two bands (p, q; or W1, W2, respectively) and is defined as

$$\begin{equation} J = \frac{1}{N} \sum _{i=1}^N \;{\rm sgn}\;(P_i) \sqrt{|{P_i}|}, \end{equation} \tag{ 2 }$$

$$\begin{equation} P_i = \delta _{p}(i)\delta _{q}(i), \end{equation} \tag{ 3 }$$

$$\begin{equation} \delta _p(i) = \sqrt{\frac{N}{N-1}} \frac{m_{p,i} - \bar{m}_{p}}{\sigma _{p,i}}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \bar{m}_{p} = \frac{1}{N}\sum _{i=1}^N m_{p,i}, \end{equation} \tag{ 5 }$$

where i is the index for each data point, N is the total number of points, sgn(P_i) is the sign of P_i, m_{p, i} is the photometric magnitude of flux measurement i in band p, and σ_{p, i} is its uncertainty. The Stetson-K index is a measure of the kurtosis of the magnitude distribution and is calculated by collapsing a single-band light curve:

$$\begin{equation} K = \frac{1}{\sqrt{N}} \frac{\sum _{i=1}^{N} |\delta (i)|}{\sqrt{\sum _{i=1}^{N} \delta (i)^2}}. \end{equation} \tag{ 6 }$$

For a pure sinusoid, K ≃ 0.9, while for a Gaussian magnitude distribution, K ≃ 0.798. This is also the scaling factor in the L-index (Equation (1)).

Our third derived feature is the magnitude ratio (MR; Kinemuchi et al. 2006). This measures the fraction of time a variable star spends above or below its median magnitude and is useful for distinguishing between variability from eclipsing binaries and pulsating variables. This is computed using the magnitude measurements m_i for an individual band and is defined as

$$\begin{equation} {\rm MR} = \frac{{\rm max}(m_i) - {\rm median}(m_i)}{{\rm max}(m_i) - {\rm min}(m_i)}. \end{equation} \tag{ 7 }$$

For example, if a variable star spends >50% of its time at near constant flux that only falls occasionally, MR ≈ 1. If its flux rises occasionally, MR ≈ 0. A star whose flux is more sinusoidal will have MR ≈ 0.5.

The remaining four features are derived from a Fourier decomposition of the W1 light curves. Fourier decomposition has been shown to be a powerful tool for variable star classification (Deb & Singh 2009; Rucinski 1993). To reduce the impact of noise and outliers in the photometric measurements, we first smooth a mid-IR light curve using a local non-parametric regression fit with Gaussian kernel of bandwidth (σ) 0.05 days. We then fit a fifth-order Fourier series to the smoothed light curve m(t), parameterized as

$$\begin{equation} m(t) = A_0 + \sum _{j=1}^5 A_j \cos [2\pi j \Phi (t) + \phi _j], \end{equation} \tag{ 8 }$$

where Φ(t) is the orbital phase at observation time t relative to some reference time t₀ and is computed using our recovered period P_r (see above) as

$$\begin{equation} \Phi (t) = \frac{t - t_0}{P_r} - \;{\rm int}\;\left(\frac{t-t_0}{P_r}\right), \end{equation} \tag{ 9 }$$

where "int" denotes the integer part of the quantity and 0 ⩽ Φ(t) ⩽ 1. The parameters that are fit in Equation (8) are the amplitudes A_j and phases ϕ_j. The quantities that we found useful for classification (see Section 4 with some guidance from Deb & Singh 2009) are the two relative phases

$$\begin{equation} \phi _{21} = \phi _2 - 2\phi _1 \end{equation} \tag{ 10 }$$

$$\begin{equation} \phi _{31} = \phi _3 - 3\phi _1, \end{equation} \tag{ 11 }$$

and the absolute values of two Fourier amplitudes: |A₂| and |A₄|.

To summarize, we have a feature vector consisting of seven metrics for each mid-IR light curve in our training sample: the recovered period P_r, the MR (Equation (7)), the Stetson-L index (Equation (1)), and the four Fourier parameters: |A₂|, |A₄|, ϕ₂₁, and ϕ₃₁. Together with available class information from the literature (the dependent variable), we constructed a data matrix consisting of 8273 labeled points ("truth" samples) in a seven-dimensional space. Section 5 describes how this data matrix is used to train and validate a machine-learned classifier.

ML methods based on classification and regression trees were popularized by Breiman (1984). Decision trees are intuitive and simple to construct. They use recursive binary partitioning of a feature space by splitting individual features at values (or decision thresholds) to create disjoint rectangular regions—the nodes in the tree. The tree-building process selects both the feature and threshold at which to perform a split by minimizing some measure of the inequality in the response between the two adjacent nodes (e.g., the fractions of objects across all known classes, irrespective of class). The splitting process is repeated recursively on each subregion until some terminal node size is reached (node size parameter below).

Classification trees are powerful non-parametric classifiers that can deal with complex non-linear structures and dependencies in the feature space. If the trees are sufficiently deep, they generally yield a small bias with respect to the true model that relates the feature space to the classification outcomes. Unfortunately, single trees do rather poorly at prediction since they lead to a high variance, e.g., as encountered when overfitting a model to noisy data. This is a consequence of the hierarchical structure of the tree: small differences in the top few nodes can produce a totally different tree and hence wildly different outcomes. Therefore, the classic bias versus variance tradeoff problem needs to be addressed. To reduce this variance, Breiman (1996) introduced the concept of bagging (bootstrap aggregation). Here, many trees (N_tree of them) are built from randomly selected (bootstrapped) subsamples of the training set and the results are then averaged. To improve the accuracy of the final averaged model, the RF method (Breiman 2001) extends the bagging concept by injecting further randomness into the tree building process. This additional randomness comes from selecting a random subset of the input features (m_try parameter below) to consider in the splitting (decision) process at each node in an individual tree. This additional randomization ensures the trees are more de-correlated prior to averaging and gives a lower variance than bagging alone, while maintaining a small bias. These details may become clearer in Section 5.4 where we outline the steps used to tune and train the RF classifier. For a more detailed overview, we refer the interested reader to Chapter 15 of Hastie et al. (2009) and Breiman & Cutler (2004).

Random forests are popular for classification-type problems and are used in many disciplines such as bioinformatics, Earth sciences, economics, genetics, and sociology. They are relatively robust against overfitting and outliers, weakly sensitive to choices of tuning parameters, can handle a large number of features, can achieve good accuracy (or minimal bias and variance), can capture complex structure in the feature space, and are relatively immune to irrelevant and redundant features. Furthermore, RFs include a mechanism to assess the relative importance of each feature in a trivial manner. This is explored in Section 5.5. For our tuning, training, and validation analyses, we use tools from the R statistical software environment.⁴ In particular, we make extensive use of the caret (Classification and Regression Training) ML package (Kuhn 2008), version 5.17-7, 2013 August.

As mentioned earlier, the RF method is relatively immune to features that are correlated with any other feature or some linear combination of them, i.e., that show some level of redundancy in "explaining" the overall feature space. However, it is recommended that features that are strongly correlated with others in the feature set be removed in the hope that a model's prediction accuracy can be ever slightly improved by reducing the variance from possible overfitting. Given that our training classes contain a relatively large number of objects (≳ 1200) and our feature set consists of only seven predictor variables, we do not expect overfitting to be an issue, particularly with RFs. Our analysis of feature collinearity is purely exploratory. Combined with our exploration of relative feature importance in Section 5.5, these analyses fall under the general topic of "feature engineering" and are considered common practice prior to training any ML classifier.

To quantify the level of redundancy or correlation among our M features, we used two methods: (1) computed the full pair-wise correlation matrix and (2) tested for general collinearity by regressing each feature on a linear combination of the remaining M − 1 features and examining the magnitude and significance of the fitted coefficients.

Figure 6 shows the pair-wise correlation matrix where elements were computed using Pearson's linear correlation coefficient ρ_ij = cov(i, j)/(σ_iσ_j) for two features i, j where cov is their sample covariance and σ_i, σ_j are their sample standard deviations. It's important to note that this only quantifies the degree of linear dependency between any two features. This is the type of dependency of interest since it can immediately allow us to identify redundant features and hence reduce the dimensionality of our problem. The features that have the largest correlation with any other feature are |A₂|, |A₄|, and L-index. Although relatively high and significant (with a <0.01% chance of being spurious), these correlations are still quite low compared to the typically recommended value of ρ ≈ 0.9 at which to consider eliminating a feature.

Zoom In Zoom Out Reset image size

The pair-wise correlation matrix provided a first crude check, although there could still be hidden collinearity in the data whereby one or more features are captured by a linear combination of the others. Our second test was therefore more general and involved treating each feature in turn as the dependent variable and testing if it could be predicted by any or all of the remaining M − 1 features (the independent variables in the regression fit). We examined the R² values (coefficients of determination) from each fit. These values quantify the proportion of the variation in the dependent variable that can be "explained" by some linear combination of all the other variables (features). The highest R² value was 0.68 and occurred when |A₂| was the dependent variable. This was not high enough to warrant removing this feature. We also explored the fitted coefficients and their significance from each linear fit. Most of them were not significant at the <5% level. We conclude that none of the features exhibit sufficiently strong correlations or linear dependencies to justify reducing the dimensionality of our feature space.

Before attempting to tune and train an RF classifier, we first partition the input training sample described in Section 3.1 into two random subsamples, containing 80% and 20% of the objects, or 6620 and 1653 objects, respectively. These are respectively referred to as the true training sample for use in tuning and training the RF model using recursive cross-validation (see below), and a test sample for performing a final validation check and assessing classification performance by comparing known to predicted outcomes. This test sample is sometimes referred to as the hold-out sample. The reason for this random 80/20 split is to ensure that our performance assessment (using the test sample) is independent of the model development and tuning process (on the true training sample). That is, we do not want to skew our performance metrics by a possibly overfitted model, however slight that may be.

In Figure 7, we compare the W1 magnitude distributions for the true training sample (referred to as simply "training sample" from here on), the test sample to support final validation, and from this, an even smaller test subsample consisting of 194 objects with W1 magnitudes ⩽9 mag. This bright subsample will be used to explore the classification performance for objects with a relatively higher S/N. The shape of the training and test-sample magnitude distributions are essentially equivalent given both were randomly drawn from the same input training sample as described above. That is, the ratio of the number of objects per magnitude bin from each sample is approximately uniform within Poisson errors. The magnitudes in Figure 7 are from the WISE All-Sky Release Catalog (Cutri et al. 2012) and derived from simultaneous point-spread function (PSF) fit photometry on the stacked single exposures covering each source from the four-band cryogenic phase of the mission only. These therefore effectively represent the time-averaged light-curve photometry. W1 saturation sets in at approximately 8 mag, although we included objects down to 7 mag after examining the quality of their PSF-fit photometry and light curves. The PSF-fit photometry was relatively immune to small amounts of saturation in the cores of sources.

Zoom In Zoom Out Reset image size

The S/N limits in Figure 7 are approximate and based on the rms noise in repeated single-exposure photometry for the non-variable source population (section IV.3.b.ii. in Cutri et al. 2012). These represent good overall proxies for the uncertainties in the light-curve measurements at a given magnitude. For example, the faintest sources in our training sample have a time-averaged W1 magnitude of ≈14.3 mag where S/N ≈14 and hence σ ≈ 1.086/14 ≈ 0.08 mag. This implies the fainter variables need to have progressively larger variability amplitudes in order to be reliably classified, e.g., with say ≳ 5σ or ≳ 0.4 mag at W1 ≳ 14 mag. This therefore sets our effective sensitivity limit for detecting and characterizing variability in the WISE single-exposure database. The paucity of faint high-amplitude variables among the known variable-source population in general explains the gradual drop-off in numbers beyond W1 ≈11.3 mag.

An overview of the RF ML algorithm was given in Section 5.1. Even though the RF method is referred to as a non-parameteric classification method, it still has a number of tuning parameters to control its flexibility. In this paper, these are (1) the number of decision trees N_tree to build from each boostrapped sample of the training set, (2) the number of features m_try to randomly select from the full set of M features to use as candidates for splitting at each tree node, and (3) the size of a terminal node in the tree, node size, represented as the minimum number of objects allowed in the final subregion where no more splitting occurs. For classification problems (as opposed to regression where the response is a multi-valued step function), Breiman (2001) recommends building each individual tree right down to its leaves where node size =1, i.e., leaving the trees "unpruned." This leaves us with N_tree and m_try. The optimal choice of these parameters depends on the complexity of the classification boundaries in the high-dimensional feature space.

We formed a grid of N_tree and m_try test values and our criterion for optimality (or figure-of-merit) was chosen to be the average classification accuracy. This is defined as the ratio of the number of correctly predicted classifications from the specific RF model to the total number of objects in all classes. This metric is also referred to as the average classification efficiency and "1 − accuracy" is the error rate. This is a suitable metric to use here since our classes contain similar numbers of objects and the overall average accuracy will not be skewed toward the outcome for any particular class. This metric would be biased for instance if one class was substantially larger than the others since it would dictate the average classification accuracy.

Fortunately, the classification accuracy is relatively insensitive to N_tree when N_tree is large (> a few hundred) and m_try is close to its optimum value. The only requirement is that N_tree be large enough to provide good averaging ("bagging") to minimize the tree-to-tree variance and bring the prediction accuracy to a stable level, but not too large as to consume unnecessary compute runtime. Therefore, we first fixed N_tree at a relatively large value of 1000, then tuned m_try using a ten-fold cross-validation on the true training sample defined in Section 5.3 and selecting the m_try value that maximized the classification accuracy (see below). Once an optimal value of m_try was found, we then explored the average classification accuracy as a function of N_tree to select an acceptable value of N_tree. Ten-fold cross-validation (or K-fold in general) entails partitioning the training sample into 10 subsamples where each subsample is labeled k = 1, 2, 3...10, then training the RF model on nine combined subsamples and predicting classifications for the remaining one. These predictions are compared to the known (true) classifications to assess classification performance. This prediction subsample is sometimes referred to as the "hold-out" or "out-of-bag" sample. Given N objects in the training sample, we iterate until every subsample k containing N/10 objects has served as the prediction data set using the model trained on the remaining T = 9N/10 objects. The final classification (or prediction) performance is then the average of all classification accuracies from all 10 iterations.

The caret package in R provides a convenient interface to train and fit an RF model using K-fold cross-validation. This calls the higher level randomForest() function, an implementation of the original Breiman & Cutler (2004) algorithm written in Fortran. We were unable to find a precise description of the algorithm implemented in tandem with cross-validation by the R caret package. Given that there is a lot happening in the training and tuning phase, we lay out the steps in Appendix A.

Figure 8 shows the average classification accuracy as a function of the trial values of m_try. The optimal value is m_try = 2 and close to the rule of thumb suggested by Breiman (2001): $m_{\rm try}\approx \sqrt{M}$, where M is the total number of features (7 here). As mentioned earlier, the classification accuracy is relatively insensitive to N_tree when N_tree is large and m_try is close to its optimum value. The results in Figure 8 assume a fixed value N_tree = 1000. Figure 9 shows the average classification accuracy as a function of N_tree for m_try = 2, 3, and 4. The achieved accuracies are indeed independent of N_tree for N_tree ≳ 400. However, to provide good tree-averaging ("bagging") and hence keep the variance in the final RF model fit as small as possible, we decided to fix N_tree at 700. This also kept the compute runtime at a manageable level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When predicting the classification for a new object with feature vector X, it is pushed down the tree. That is, it is assigned the label of the training sample in the terminal node it ends up in. This procedure is iterated over all N_tree trees in the ensemble, and the mode (or majority) vote of all trees is reported as the predicted class. However, instead of the winning class, one may want to quote the probabilities that an object belongs to each respective class. This allows one to make a more informed decision. The probability that a new object with feature vector X belongs to some class C_j where j = 1, 2, 3, ... is given by

$$\begin{equation} P(\mathbf {X}|C_j) = \frac{1}{N_{\rm tree}} \sum _{i=1}^{N_{\rm tree}}I(p_i=C_j|\mathbf {X}), \end{equation} \tag{ 12 }$$

where I(p_i = C_j|X) is an indicator function defined to be 1 if tree i predicts class C_j and 0 otherwise. In other words, for the RF classifier, the class probability is simply the fraction of N_tree trees that predicted that class. It is important to note that the probability computed via Equation (12) is a conditional class probability and only has meaning relative to the probabilities of obtaining the same features X conditioned on the other contending classes. That is, we say an object with features X is relatively more likely to have been generated by the population of objects defined by class C₁ than classes C₂ or C₃, etc.

The definition in Equation (12) should not be confused with the posterior probability that the object belongs to class C_j in an "absolute" sense, i.e., as may be inferred using a Bayesian approach. This can be done by assuming some prior probability P(C_j) derived for example from the proportion that each class contributes to the observable population of variable sources. The probability in this case would be:

$$\begin{equation} P(C_j|\mathbf {X}) = \frac{P(\mathbf {X}|C_j) P(C_j)}{P(\mathbf {X})}, \end{equation} \tag{ 13 }$$

where P(X) is the normalization factor that ensures the integral of P(C_j|X) over all C_j is 1 and P(X|C_j), the "likelihood," is given by Equation (12). Unfortunately, plausible values for the priors P(C_j) are difficult to derive at this time since the relative number of objects across classes in our training sample are likely to be subject to heavy selection biases, e.g., brought about by both the WISE observational constraints and heterogeneity of the input optical variability catalogs used for the cross-matching. The current relative proportions of objects will not represent the true mix of variable sources one would observe in a controlled flux-limited sample according to the WISE selection criteria and sensitivity to each class. This Bayesian approach will be considered in future as classification statistics and selection effects are better understood. For our initial classifications, we will quote the relative (conditional) class probabilities defined by Equation (12) (see Section 6 for details).

An initial qualitative assessment of the separability of our three broad variability classes according to the seven mid-IR light-curve features was explored in Section 4. Given the relatively high dimensionality of our feature space, this separability is difficult to ascertain by looking at pairwise relationships alone. Our goal is to explore class separability using the full feature space in more detail. This also allows us to identify those features that best discriminate each class as well as those that carry no significant information, both overall and on a per-class basis.

RFs provide a powerful mechanism to measure the predictive strength of each feature for each class, referred generically to as feature importance. This quantifies, in a relative sense, the impact on the classification accuracy from randomly permuting a feature's values, or equivalently, forcing a feature to provide no information, rather than what it may provide on input. This metric allows one to determine which features work best at distinguishing between classes and those that physically define each class. Feature importance metrics can be additionally generated in the RF training/tuning phase using the "hold-out" (or "out-of-bag") subsamples during the cross-validation process (Section 5.4), i.e., excluded from the bootstrapped training samples. The prediction accuracies from the non-permuted-feature and permuted-feature runs are differenced then averaged over all the N_tree trees and normalized by their standard error ($\hat{\sigma }/\sqrt{N_{\rm tree}}$). The importance metrics are then placed on a scale of 0–100 where the "most important" metric is assigned a value of 100 for the class where it is a maximum. The intent here is that features that lead to large differences in the classification accuracy for a specific class when their values are randomly permuted are also likely to be more important for that class. It is worth noting that even though this metric is very good at finding the most important features, it can give misleading results for highly correlated features which one might think are important. A feature could be assigned a low RF importance score (i.e., with little change in the cross-validation accuracy after its values are permuted) simply because other features that strongly correlate with it will stand in as "surrogates" and carry its predictive power.

Figure 10 illustrates the relative importance of each feature for each class using the predictions from our initially tuned RF classifier. Period is the most important feature for all classes. This is no surprise since the three classes are observed to occupy almost distinct ranges in their recovered periods (Figure 1), therefore providing good discriminative and predictive power. In general, not all features are equally as important across classes. For example, the relative phase ϕ₂₁ and L-index are relatively weak predictors on their own for the W UMa and Algol eclipsing binaries respectively while they are relatively strong predictors for RR Lyr pulsating variables. In practice, one would eliminate the useless features (that carry no information) across all classes and retrain the classifier. However, since all features have significant predictive power for at least one class in this study, we decided to retain all features.

Zoom In Zoom Out Reset image size

We explored the impact on the per-class prediction accuracy of randomly permuting the ϕ₂₁ and L-index feature values, i.e., so they provide no useful information. The RF classifier was retrained and the final validation test sample defined in Section 5.3 was used to assess the classification accuracy. This provided more of an "absolute" measure of the importance of these features than in Figure 10. We found that if either ϕ₂₁ or L-index were forced to be useless, the classification accuracies for the Algol and W UMa classes were not significantly altered. However, the accuracy dropped by ≈4.1% and ≈3.1% for the RR Lyr class by forcing these features to be useless respectively. If both ϕ₂₁ and L-index were forced to be useless, the change in classification accuracies for the Algol and W UMa classes were still insignificant (dropping by ≲ 1.5%), but the drop for the RR Lyr class was ≈7.7%. This not only confirms the results in Figure 10, but also the fact that RF classifiers are relatively immune to features that carry little or no class information.

We validate the overall accuracy of the RF classifier that was fit to the training sample by predicting classifications for the two test samples defined in Section 5.3 (Figure 7) and comparing these to their "true" (known) classifications. These test samples are independent of the training sample and hence allow an unbiased assessment of classifier performance. This was explored by computing the confusion matrix across all classes. The confusion matrix for our largest test sample (consisting of 1653 objects to W1 ∼14 mag) is shown in Figure 11. The quantities therein represent the proportion of objects in each true (known) class that were predicted to belong to each respective class, including itself. The columns are normalized to add to unity. When compared to itself (i.e., a quantity along the diagonal going from top-left to bottom-right in Figure 11), it is referred to as the sensitivity in ML parlance. It is also loosely referred to as the per-class classification accuracy, efficiency (e.g., as in Section 5.5), or true positive rate. We obtain classification efficiencies of 80.7%, 82.7%, and 84.5% for Algols, RR Lyrae, and W UMa type variables respectively. The overall classification efficiency, defined as the proportion of all objects that were correctly predicted (irrespective of class) is ≈82.5%. The corresponding 95% confidence interval (from bootstrapping) is 80.5%–84.3%, or approximately ±2% across all three classes.

Zoom In Zoom Out Reset image size

For comparison, Richards et al. (2011) obtained an overall classification efficiency of ≈77.2% on a 25 class data set of 1542 variable stars from the OGLE and Hipparcos surveys. However, if we isolate their Algol (predominately β Lyrae), RR Lyrae (all types), and W UMa statistics, we infer an overall classification efficiency of ≈88.4%, implying an improvement of ≈6% over our estimate for WISE variables. This difference is likely due to their higher quality, longer time span optical light curves—specially selected to have been well studied in the first place. Nonetheless, our classification performance is still good given the WISE cadence, sparsity and time span of observations, and possible uncertainties in classifications from the literature used to validate the predictions.

The off-diagonal quantities of the confusion matrix in Figure 11 can be used to compute the reliability (or equivalently the purity, specificity, or "1 − false positive rate") for a specific class. That is, the proportion of objects in all other classes that are correctly predicted to not contaminate the class of interest. This can be understood by noting that the only source of unreliability (or contamination) in each class are objects from other classes. For example, the false positive rate (FPR) for the Algol class is

$$\begin{eqnarray*} &&{{\displaystyle {\rm FPR} = \frac{0.112\times 347 + 0.118\times 587}{347\times \left(0.112+0.827+0.061\right)+ 587\times \left(0.118+0.037+0.845\right)}\approx 0.116 }}\\ \end{eqnarray*}$$

and hence its purity is 1 − FPR ≈ 88.4%. Similarly, the purity levels for the RR Lyrae and W UMa classes from Figure 11 are 96.2% and 87.6% respectively. For comparison, Richards et al. (2011) obtain purity levels of up to 95% for these and most other classes in their study.

For the smaller, higher S/N test sample of 194 objects with W1 magnitudes ⩽9 (Figure 7), the classification accuracy for the Algol class improves to ≈89.7%, compared to 80.7% for the large test sample. However for the RR Lyr class, the classification accuracy drops to 55.5% (from 82.7%) and for the W UMa class, it drops to 79.3% (from 84.5%). In general, we would have expected an increase in classification accuracy across all classes when only higher S/N measurements, and hence objects with more accurately determined features are used. This indeed is true for the Algol class which appears to be the most populous in this subsample with 127 objects. The drop in classification performance for the other two classes can be understood by low number statistics with only 9 RR Lyr and 58 W UMa objects contributing. Their sampling of the feature space density distribution in the training set for their class is simply too sparse to enable reliable classification metrics to be computed using ensemble statistics on the predicted outcomes. In other words, there is no guarantee that most of the nine RR Lyr in this high S/N test sample would fall in the densest regions of the RR Lyr training model feature space so that they can be assigned high enough RF probabilities to be classified as RR Lyr.

The primary output from the RF classifier when predicting the outcome for an object with given features is a vector of conditional class likelihoods as defined by Equation (12). By default, the "winning" class is that with the highest likelihood. A more reliable and secure scheme to assign the winning class will be described in Section 6. Distributions of all the classification probabilities for our largest test sample are shown in Figure 12. These probabilities are conditioned on the winning class that was assigned by the RF classifier so that histograms at the high end of the probability range in each class-specific panel correspond to objects in that winning class. The spread in winning class probabilities is similar across the three classes, although the Algol class has slightly more mass at P ≳ 0.7 (Figure 12(a)). This indicates that the seven-dimensional feature space sample density is more concentrated (or localized) for this class than for the other classes.

Zoom In Zoom Out Reset image size

Figure 13 shows the receiver operating characteristic (ROC) curves for each class in our largest test sample. These are generated by thresholding on the classification probabilities of objects in each class (i.e., with P > 0, P > 0.02, P > 0.04, ..., P > 0.98 from left to right in Figure 13), then computing the confusion matrix for each thresholded subclass. The true positive rate (TPR or classification accuracy) and false positive rate (FPR or impurity) were then extracted to create the ROC curve. The trends for these curves are as expected. Given the class probability quantifies the degree of confidence that an object belongs to that class, the larger number of objects sampled to a lower cumulative probability level will reduce both the overall TPR and FPR. That is, a smaller fraction of the truth is recovered, but the number of contaminating objects (false positives) from other classes does not increase much and the larger number of objects in general keeps the FPR relatively low. The situation reverses when only objects with a higher classification probability are considered. In this case there are fewer objects in total and most of them agree with the true class (higher TPR). However, the number of contaminants is not significantly lower (or does not decrease in proportion to the reduced number of objects) and hence the FPR is slightly higher overall.

Zoom In Zoom Out Reset image size

It is also interesting to note that even though the full test sample confusion matrix in Figure 11 indicates that W UMa objects have the highest classification accuracy (at 84.5%—corresponding to far left on the ROC curve in Figure 13), this is overtaken by RR Lyrae at internal probability thresholds of P ≳ 0.1 where the classification accuracy (TPR) becomes >86%. This however is at the expense of an increase in the FPR to >12%. Therefore, the ROC curves contain useful information for selecting (class-dependent) classification probability thresholds such that specific requirements on the TPR and FPR can be met.

We compare the performance of the RF classifier trained above to other popular machine-learned classifiers. The motive is to provide a cross-check using the same training data and validation test samples. We explored artificial neural networks (NNET), kNN, and SVM. A description of these methods can be found in Hastie et al. (2009). The R caret package contains convenient interfaces and functions to train, tune, test, and compare these methods (Kuhn 2008). More ML methods are available, but we found these four to be the simplest to set up and tune for our problem at hand. Furthermore, these methods use very dissimilar algorithms and thus provide a good comparison set. Parameter tuning was first performed automatically using large grids of test parameters in a ten-fold cross-validation (defined in Section 5.4); then the parameter ranges were narrowed down to their optimal ranges for each method using grid sizes that made the effective number of computations in training approximately equal across methods. This enabled a fair comparison in training runtimes.

The classification accuracies (or efficiencies), runtimes, and the significance of the difference in mean accuracy relative to the RF method are compared in Table 1. The latter is in terms of the p-value of obtaining an observed difference of zero (the null hypothesis) by chance according to a paired t-test. It appears that the NNET method performs just as well as the RF method in terms of classification accuracy, although the RF method has a slight edge above the others. This can be seen in Figure 14 where the overall distributions in accuracy are compared. Aside from the similarity in classification performance between NNET and RF, the added benefits of the RF method, e.g., robustness to outliers, flexibility and ability to capture complex structure, interpretability, relative immunity to irrelevant and redundant information, and simple algorithms to measure feature importance and proximity for supporting active learning frameworks (Section 6.1), makes RF our method of choice.

Zoom In Zoom Out Reset image size

It is known that RF classifiers trained using supervised methods perform poorly outside their "learned boundaries," i.e., when extrapolating beyond their trained feature space. The RF training model constructed in Section 5.4 was tuned to predict the classifications of only three broad classes: Algol, RR Lyrae, and W UMa—the most abundant types that could be reliably identified given the WISE sensitivity and observing cadence. Furthermore, this model is based on confirmed variables and classifications from previous optical surveys (Section 3.1) which no doubt contain some incorrect labels, particularly since most of these studies also used some automated classification scheme. Therefore, our initial training model is likely to suffer from sample selection bias whereby it will not fully represent all the variable types that WISE can recover or discover down to fainter flux levels and lower S/Ns (Figure 7). Setting aside the three broad classes, our initial training model will lead to biased (incorrect) predictions for other rare types of variables that are close to or distant from the "fuzzy" classification boundaries of the input model.

Figure 15 illustrates some of these challenges. Here we show example W1 and W2 light curves for a collection of known variables from the literature (including one "unknown") and their predicted class probabilities using our input training model. The known short-period Cepheid (top left) would have its period recovered with good accuracy given the available number of WISE exposures that cover it. However, it would be classified as an Algol with a relatively high probability. That's because our training sample did not include short-period Cepheids. Its period of ∼5.6 days is at the high end of our fitted range (Figure 1) and overlaps with the Algol class. For the given optimum observation time span covered by WISE, the number of short-period Cepheids after cross-matching was too low to warrant including this class for reliable classification in future. Better statistics at higher ecliptic latitudes (where observation time spans are longer) are obviously needed. The known Algol and one of the two known RR Lyrae in Figure 15 are securely classified, although the known W UMa achieves a classification probability of only 0.535 according to our training model. This would be tagged as "unknown" if the probability threshold was 0.6 for instance. The two lower S/N objects on the bottom row (a known RR Lyra and a fainter variable we identify as a possible RR Lyra "by eye") would also be classified as "unknown" according to our initial model, even though their light curves can be visually identified as RR Lyrae. This implies that when the S/N is low and/or the result from an automated classification scheme is ambiguous (following any refinement of the training model; see below), visual inspection can be extremely useful to aid the classification process.

Zoom In Zoom Out Reset image size

We need a mechanism that fills in the undersampled regions of our training model but also improves classification accuracies for the existing classes. Richards et al. (2012, and references therein) presented methods to alleviate training-sample selection biases and we use their concepts (albeit slightly modified) to optimize and extend classifications for the WVSDB. These methods fall under the general paradigm of semi-supervised or active learning whereby predictions and/or contextual follow-up information for new data is used to update (usually iteratively) a supervised learner to enable more accurate predictions. Richards et al. (2012) were more concerned with the general problem of minimizing the bias and variance of classifiers trained on one survey for use on predicting the outcomes for another. Our training sample biases are not expected to be as severe since our training model was constructed more-or-less from the same distribution of objects with properties that we expect to classify in the long run. Our goal is simply to strengthen predictions for the existing (abundant) classes as well as explore whether more classes can be teased out as the statistics improve and more of the feature space is mapped. Our classification process will involve at least two passes, where the second pass (or subsequent passes) will use active learning concepts to refine the training model. A review of this classification process is given in Appendix B.