# Photometric Classifications of Evolved Massive Stars:
Preparing for the Era of Webb and Roman with Machine Learning

Trevor Z. Dorn-Wallenstein University of Washington Astronomy Department
Physics and Astronomy Building, 3910 15th Ave NE
Seattle, WA 98105, USA  James R. A. Davenport University of Washington
Astronomy Department
Physics and Astronomy Building, 3910 15th Ave NE
Seattle, WA 98105, USA  Daniela Huppenkothen SRON Netherlands Institute for
Space Research,
Sorbonnelaan 3, 3584 CA Utrecht
The Netherlands  DIRAC Institute,
Department of Astronomy,
University of Washington, 3910 15th Ave NE
Seattle, WA 98105, USA  The University of Washington eScience Institute,
The Washington Research Foundation Data Science Studio,
University of Washington,
Seattle, WA 98105, USA  Emily M. Levesque University of Washington Astronomy
Department
Physics and Astronomy Building, 3910 15th Ave NE
Seattle, WA 98105, USA

###### Abstract

In the coming years, next-generation space-based infrared observatories will
significantly increase our samples of rare massive stars, representing a
tremendous opportunity to leverage modern statistical tools and methods to
test massive stellar evolution in entirely new environments. Such work is only
possible if the observed objects can be reliably classified. Spectroscopic
observations are infeasible with more distant targets, and so we wish to
determine whether machine learning methods can classify massive stars using
broadband infrared photometry. We find that a Support Vector Machine
classifier is capable of coarsely classifying massive stars with labels
corresponding to hot, cool, and emission line stars with high accuracy, while
rejecting contaminating low mass giants. Remarkably, 76% of emission line
stars can be recovered without the need for narrowband or spectroscopic
observations. We classify a sample of ${\sim}2500$ objects with no existing
labels, and identify fourteen candidate emission line objects. Unfortunately,
despite the high precision of the photometry in our sample, the heterogeneous
origins of the labels for the stars in our sample severely inhibits our
classifier from distinguishing classes of stars with more granularity.
Ultimately, no large and homogeneously labeled sample of massive stars
currently exists. Without significant efforts robustly classify evolved
massive stars — which is feasible given existing data from large all-sky
spectroscopic surveys — shortcomings in the labeling of existing data sets
will hinder efforts to leverage the next-generation of space observatories.

††journal: ApJ††software: Astropy v3.2.2 (Astropy Collaboration et al., 2013;
The Astropy Collaboration et al., 2018), Astroquery v0.3.10 (Ginsburg et al.,
2019), Matplotlib v3.1.1 (Hunter, 2007), makecite (Price-Whelan et al., 2018),
NumPy v1.17.2 (Van Der Walt et al., 2011), Pandas v0.25.1 (McKinney, 2010),
Python 3.7.4, Scikit-learn v0.21.3 (Pedregosa et al., 2011), Scipy v1.3.1
(Jones et al., 2001)

##  1 Introduction

Evolved massive stars are observed in a menagerie of exotic evolutionary
phases. While the challenge of connecting these states with a self-consistent
theory of stellar evolution has seen rapid advancement since the original
introduction of the “Conti Scenario” (Conti et al., 1983), the effects of
rotation, magnetic fields, internal mixing processes, and binary interactions
on the evolution of massive stars are still the subject of much theoretical
effort (e.g., Ekström et al., 2012; Eldridge et al., 2017). While individual
massive stars can be used as precision probes of these processes, ensembles of
evolved massive stars can also significantly constrain stellar evolution. This
can be done by comparing the integrated spectra of massive stars (e.g.
Levesque et al., 2012), or by studying the detailed makeup of resolved
populations of massive stars (Dorn-Wallenstein & Levesque, 2018, 2020; Stanway
et al., 2020).

The latter approach requires large and accurately-classified samples of
evolved massive stars that will be achievable in the coming years with the
launch of the James Webb Space Telescope (Webb) and the Nancy Grace Roman
Space Telescope (Roman). Among the instrumentation on Webb and the proposed
instrumentation for Roman are photometers equipped with filters spanning a
broad wavelength baseline from 0.5 to 28 $\mu$m. The resolution of Webb will
allow us to identify and study in detail individual luminous stars to great
distances (e.g. Jones et al., 2017), while the impressive 0.218 deg2 field of
view of Roman will allow us to efficiently survey nearby galaxies in a small
number of pointings (Spergel et al., 2013). Combined, observations from both
missions will give astronomers access to precise infrared measurements of vast
numbers of evolved massive stars. But without sophisticated methods of
identifying and classifying these stars, the science return afforded by such a
large increase in expected sample sizes will be significantly reduced.

Classification of stars from broadband photometry is often done by adopting
simple linear cuts in color-magnitude space (e.g., Massey et al., 2006, 2009),
and — most critically — do not include rare emission line objects, whose
classification requires dedicated narrow-band surveys (sometimes with custom-
designed filters, e.g. Neugent et al. 2018b), often accompanied by follow-up
spectroscopy, both of which require extensive telescope time. These objects
are often the post-main sequence evolved states of massive stars, in which the
effects of rotation, binary interactions, and chemical mixing are the most
pronounced; the stars that place the most valuable constraints on unknown
stellar physics are also the hardest to detect via traditional means.
Therefore, it is worthwhile to determine whether there are alternative ways to
classify massive stars that avoid using traditional and expensive methods.

At present, we can mimic the observing capabilities of Webb and Roman by
combining data from Gaia (which has a red-optical bandpass), the Two Micron
All Sky Survey (2MASS, Skrutskie et al. 2006, near-infrared), and the Wide
Field Infrared Survey Explorer (WISE, Wright et al. 2010, mid-infrared). WISE
provides the additional benefit of having scanned the sky approximately every
six months, yielding lightcurves spanning a $\sim 7$-year baseline from which
we can extract variability metrics for most stars observed. While Roman and
Webb will not be observing the entire sky in this fashion, determining whether
variability can aid in the classification of evolved massive stars will
determine whether observers should seek repeated observations of a stellar
population.

We wish to determine whether we can

  1. 1.

Assemble a sample of evolved massive stars with available classifications as a
training data set,

  2. 2.

Construct a machine learning classifier that can reject low mass red
contaminants and identify likely emission line objects in order to optimise
available telescope time on the most promising targets,

  3. 3.

Determine whether variability metrics estimated from WISE lightcurves can aid
in these tasks, and

  4. 4.

Determine which photometric bandpasses and variability metrics contribute the
most to making accurate classifications.

Here we utilize a support vector machine classifier (SVC) trained only on
broad-band photometry and simple metrics derived from WISE lightcurves to
classify a large sample of evolved massive stars. We describe our sample
selection and labeling method in §2. §3.1 details the calculation of the
simple metrics derived from the WISE lightcurves, and describes the overall
behavior of the stars in our sample. We explain our classification algorithm,
discuss its successes and shortcomings in §4, and apply it to a training
sample of 2500 stars before presenting our recommendations and concluding in
§5.

##  2 Sample Selection & Labeling

For any machine learning algorithm, a high-quality training set with accurate
labels is necessary. The second data release (DR2) of the Gaia mission (Gaia
Collaboration et al., 2018a) contains precise photometry in three bands ($G$,
$G_{BP}$, and $G_{RP}$) and geometric parallaxes ($\varpi$) for 1.3 billion
stars in the Galaxy and Magellanic Clouds. Because the parallax measurements
suffer from some systematics (Lindegren et al., 2018), and many objects have
high fractional errors ($\sigma_{\varpi}/\varpi$) or negative measured
parallax, Bailer-Jones et al. (2018) calculated Bayesian distance estimates
for the majority of stars in Gaia DR2, using a prior based on the spatial
distribution of stars in the Galaxy. Figure 1 shows the difference between the
distance inferred by Bailer-Jones et al. (2018), $r_{est}$, and a naive
distance derived by inverting the reported Gaia measurements of $\varpi$ for
$\sim 10,000$ putative massive stars (as described below). The dashed line
indicates where $r_{est}=1/\varpi$. While the two distance estimates are
roughly consistent for nearby stars, more distant stars are biased much
further away in the naive distance estimates.

We first perform a cross-match between the Bailer-Jones et al. (2018) catalog
and the existing cross-match between Gaia DR2 and the ALLWISE data release.
ALLWISE (Cutri et al., 2013) contains photometry in four mid-infrared (MIR)
bands — $W1$ (3.4 $\mu$m), $W2$ (4.6 $\mu$m), $W3$ (12 $\mu$m), and $W4$ (22
$\mu$m) — derived from co-added images obtained during the original WISE
mission, as well as $W1$ and $W2$ images obtained in the post-cryogenic
NEOWISE mission (Mainzer et al., 2011). We select all stars with successful
distance estimates (i.e., where ${\tt result\\_flag}=1$ if the distance
estimate is the mode of the posterior distribution, 2 if it is the median, and
0 for a failed estimate, see Bailer-Jones et al. 2018 for more details) that
satisfy

$\displaystyle\begin{split}M_{G}&=G-5\log{r_{est}}+5\leq-1.5,\\\
W1&<14.\end{split}$ (1)

Since Bailer-Jones et al. (2018) used a Galactic prior, stars in the
Magellanic Clouds have distances that are considerably underestimated. Thus,
we also match the catalog in Gaia Collaboration et al. (2018b) to the
ALLWISE/Gaia cross-match, and select stars with $W1<14$ and $M_{G}\leq-1.5$,
assuming distance moduli of 19.05/18.52 for the SMC/LMC respectively (Kovács,
2000a, b) and combining the two cross-matches while dropping duplicate stars.
This results in a total of 452,283 stars.

Figure 1: Distance from Bailer-Jones et al. (2018) versus distance inferred
via inverting the reported $\varpi$ from Gaia DR2 for $\sim 10000$ putative
massive stars. The dashed line shows where $r_{est}=1/\varpi$.

We then estimate the reddening in the Gaia bandpasses using the published
estimate for $A_{G}$ from Gaia DR2, and coefficients from (Malhan et al.,
2018) to calculate $E(G_{BP}-G_{RP})$. For Galactic stars without $A_{G}$
estimates, we assume $A_{G}=0$, and for stars in the Magellanic Clouds, we
assume the average value of $A_{G}$ and $E(G_{BP}-G_{RP})$ using $R_{V}$
measurements from Gordon et al. (2003) and $E(B-V)$ from Massey et al. (2007).
Using these quanitites, we calculate the intrinsic $G_{BP}-G_{RP}$ and $M_{G}$
for all stars.

Figure 2: Left: Density of stars selected from the Gaia database on the sky.
Intensity of the colormap corresponds to the logarithm of the number of stars
in each bin. Right: Gaia CMD for stars selected from the Gaia DR2 database
brighter than $M_{G}=-2.5$. Galactic stars are in blue, LMC stars are in
orange, and SMC stars are in green. The solid, dashed, and dotted lines
represent our minimum-luminosity criteria to select massive stars in the
Galaxy, LMC, and SMC respectively.

We can then construct color-magnitude diagrams (CMDs) in the Gaia filters,
which we can use to select massive stars — i.e., stars with initial mass
$M_{i}\geq 8M_{\odot}$. We use the MESA Isochrones & Stellar Tracks (MIST,
Dotter 2016; Choi et al. 2016; Paxton et al. 2011, 2013, 2015) isochrones with
metallicity $[Fe/H]=0,-0.5,-1$ for the Galaxy, LMC, and SMC respectively, and
rotation speed relative to critical of $v/v_{crit}=0.4$. We then selected the
faintest isochrone point of any age between $10^{5}$ and $10^{10.3}$ yr in
steps of 0.05 dex with $M_{i}\geq 8$ M⊙ in 100 equally-spaced bins in the
range $-0.25\leq G_{BP}-G_{RP}\leq 3$. These isochrone points form a boundary
in the Gaia CMD that represents the faintest luminosities reached by any
massive star at any point during its evolution, and no fainter massive stars
are expected to be found — note that many isochrone points with $M_{i}<8$ M⊙
lie above this boundary, so our sample is not constructed to be free of
contamination. The left panel of Figure 2 shows the logarithmic density on the
sky of all stars selected from the Gaia DR2 database; the Galactic plane and
Magellanic Clouds are clearly visible. The right panel shows the portion of
the Gaia CMD brighter than $M_{G}=-2.75$, where blue (orange, green) points
are individual stars in the Galaxy (LMC, SMC). The solid (dashed, dotted)
black line show the MIST luminosity threshold for the Galaxy (LMC, SMC). Note
that the thresholds accurately capture the slope of the main sequence for all
three galaxies, as well as the $G_{BP}-G_{RP}$ color corresponding to the
Hayashi limit.

We select all stars brighter than the corresponding luminosity threshold for
their host galaxy, resulting in 9784 objects. From this sample, we select all
stars fainter than the saturation limit in $W1$ (8) and $W2$ (7) with valid
measurements listed in the ALLWISE catalog for the three bluest WISE bands
(excluding $W4$, where the signal-to-noise is often poor). We also convert the
$W1$ and $W2$ magnitudes (and uncertainties) to fluxes, and filter for stars
with signal-to-noise ratio greater than 3. This results in a final sample of
6484 objects. Table 1 lists the names, coordinates, host galaxies, distances
from Bailer-Jones et al. (2018), and Gaia photometry for these stars. We query
${\tt Vizier}$ (Ochsenbein et al., 2000) using ${\tt astroquery}$ to download
$JHK_{s}$ photometry from the 2-micron All Sky Survey (2MASS, Skrutskie et al.
2006) for all stars. We also query SIMBAD (Wenger et al., 2000) and download
the common name (${\tt MAIN\\_ID}$), spectral type (contained in the ${\tt
MK\\_Spectral\\_Type}$ and ${\tt SP\\_Type}$ fields), and object type (${\tt
OType}$) for each star, the latter two of which we use to assign labels.

Table 1: Common names, coordinates, host galaxies, and Gaia measurements of
6,484 putative massive stars, ordered by Right Ascension. $r_{est}$ from
Bailer-Jones et al. (2018) is given for Galactic stars. Listed values of $G$
and $G_{BP}-G_{RP}$ are uncorrected for extinction.

Common Name | R.A. [deg] | Dec [deg] | Host Galaxy |  $r_{est}$ [kpc] |  $G$ [mag] |  $A_{G}$ [mag] |  $G_{BP}-G_{RP}$ [mag]
---|---|---|---|---|---|---|---
HD 236270 | 0.17442287 | 55.72245665 | Milky Way | 2.162 | 9.073082 | 0.942 | 0.24057600
LS I +64 10 | 0.38838658 | 64.51219232 | Milky Way | 5.305 | 11.545427 | 1.3305 | 0.56807750
LS I +60 69 | 0.55506135 | 60.43828347 | Milky Way | 5.866 | 11.854865 | 1.234 | 0.61282300
BD+62 2353 | 0.59458742 | 62.90087875 | Milky Way | 5.243 | 9.809305 | 0.485 | 0.37020600
HD 73 | 1.40408512 | 43.40139506 | Milky Way | 1.869 | 8.1893425 | 0.1245 | -0.14992300
HD 240496 | 1.42175475 | 58.49541068 | Milky Way | 2.499 | 9.698675 | 1.5533 | 0.68225500
WISE J000559.28-790653.3 | 1.49713706 | -79.11483482 | SMC | - | 13.9371195 | 0.21266685 | 1.04591900
LS I +59 30 | 1.70503555 | 59.85955733 | Milky Way | 4.006 | 10.856779 | 1.188 | 0.50059550
BD+57 2870 | 1.82960982 | 58.33785301 | Milky Way | 3.893 | 9.841114 | 1.4205 | 0.82356400
BD+62 1 | 1.88805102 | 63.08030731 | Milky Way | 2.893 | 10.293984 | 1.2973 | 0.52861900

Note. — This table is published in its entirety in the machine-readable
format. A portion is shown here for guidance regarding its form and content.

###  2.1 Label Assignment

Figure 3: The makeup of our sample of massive stars. Note that the sample is
dominated by OBA stars and cool supergiants, and non-OBAe emission line stars
are the rarest stars in our sample. For readability, we have used a
logarithmic y-axis to display our sample statistics. Note that in practice,
differences in the number of stars per class are much larger than they might
appear here.

For our final sample of $\sim 6500$ stars, we wish to assign the best
available estimate of its evolutionary state. These labels can be used to
compare to the predictions of stellar population models. Note that these
evolutionary states (which are theoretical concepts) are mostly tied to
spectral appearance (which is an observable quantity). Therefore we are
assuming that, e.g., all stars with Wolf-Rayet spectra are in the same
evolutionary state (namely, the descendants of massive stars with high
luminosities that have lost their envelopes via strong winds), and that all
stars in that evolutionary state are observed as Wolf-Rayet stars. We know
that at least the former isn’t true, as some stars lose their envelopes due to
interactions with a binary companion (Eldridge et al., 2017), and the latter
is also questionable since such stars may or may not appear similar to
classical Wolf-Rayets (Götberg et al., 2018). Nevertheless, we assume that the
assigned labels are a reasonable approximation for a star’s evolutionary state
with this caveat in mind.

At the present, a database of homogeneously classified massive stars does not
currently exist. While all-sky spectroscopic surveys have observed many
massive stars, the machine learning pipelines that produce the effective
temperatures, surface gravities, and chemical compositions that would allow us
to accurately classify our sample do not cover the parameter regime in which
massive stars reside (e.g. García Pérez et al., 2016). As a result, we must
use the heterogeneous classification data available on SIMBAD. For each star,
we apply a decision tree that results in the star receiving one unique label.
Figure 4 shows a flowchart that summarizes our labeling scheme. Note that this
process highly tailored to this dataset, and some branches in the decision
tree serve only to accurately label very small numbers of stars with unique
spectral types (e.g., spectroscopically peculiar stars or X-ray binaries).
Deriving labels for known massive stars using existing sources is not trivial,
and our labelling scheme would be entirely different if a large sample of
massive stars with well-measured temperatures, surface gravities, and chemical
abundances was available.

We first use the common name and Gaia ${\tt source\\_id}$ of the star to
determine if the star belongs to the catalog of confirmed Luminous Blue
Variables (LBVs) presented in Richardson & Mehner (2018). Non-LBVs are
classifed as WR stars if “W” is in the spectral type, or the SIMBAD ${\tt
OType}$ is “*WR”. Non-WRs with “K” or “M” in their spectral type are
classified as either Red Supergiants (RSGs) or “C/S/Giant” if their SIMBAD
${\tt SP\\_Type}$ contains “III” — we keep all such low-mass contaminants in
our sample as distinguishing between RSGs and luminous low-mass giants is
still a difficult problem (Massey et al., 2009; Yang et al., 2019; Neugent et
al., 2020). The resulting sample of RSGs is pure; of the five RSGs that don’t
have luminosity class I, all of them are luminosity class Ia-II, Ib-II, or
Iab-II, which are consistent with bona fide RSGs (Levesque et al., 2005). This
is a good test of the Gaia DR2 parallaxes and Bailer-Jones et al. (2018)
distances, as cool subgiants and dwarfs with luminosity class IV or V would
have been erroneously classified as RSGs using our criteria had they been
included in our sample due to bad distance and $M_{G}$ measurements.

Non-RSGs with “F” or “G” in their spectral type are classified as Yellow
Supergiants (YSGs). However, this includes a number of blue stars. Further
inspection of these stars reveals a number of objects whose MK_Spectral_Type
field indicates these are hot stars, with spectral type O, B, or A, which we
classify as such (see below). Eight low-mass yellow stars are also included in
our sample. Stars with “III” or “V” in their spectral types are classified as
Yellow Dwarfs. We note that for extragalactic samples, foreground dwarfs can
usually be filtered based on proper motions, while dwarfs belonging to the
stellar population under study can be excluded just based on their apparent
magnitudes. Nonetheless, we retain this class to avoid confusion between these
stars and true YSGs. All YSGs that aren’t hot stars or dwarfs keep their YSG
label.

Of the objects that have not yet been classified, stars with “[e]” in their
spectral type are classified as OB[e] stars, while non-OB[e] stars with
spectral types including an “e” (without brackets, and without “pec” in their
spectral type) are classified as OBAe stars. If the star is not yet classified
and O, B, or A are in the spectral type with no additional information, they
are classified as generic OBA stars. OBA stars with “III” or “IV” in their
spectral type are classified as evolved OBA stars, stars with “V” in their
spectral type are classified as OBA main sequence stars, and finally, stars
with “I” in their spectral type are labeled as OBA supergiants. All stars that
have not been assigned a label at this stage are either C/S stars (which are
assigned the C/S/Giant class), stars labeled only as Variables in SIMBAD
(e.g., LPVs, semi-regular variables, or just variables, without other spectral
information) which are assigned the Miscellaneous Variable classification, or
stars with no identifying information/no confirmed designation (e.g., the
SIMBAD ${\tt OType}$ is “Star” or contains “Candidate”) which are classified
as Unknown/Candidate.

Figure 4: Flowchart illustrating the process by which stars are assigned
labels, as described in text. Each star begins in the top left and is assigned
a label by following a series of binary decisions. This process is complex,
and demonstrates the difficulty in deriving useful labels for massive stars.

Finally, we include an “Is Binary” flag for all stars, which is 1 for stars
classified as Eclipsing or Spectroscopic Binaries, High Mass X-ray Binaries,
or Ellipsoidal Variables, or if they have a compound spectral type (e.g., WN8
+ O6V),111This does not include stars with the “OB+” spectral type, which is
an outdated class that describes OB stars with weaker absorption lines that
would now be classified as OB supergiants. and 0 otherwise; 102 stars are
flagged as binaries. Because photometry of binary systems can be misleading
(Neugent et al., 2018a), and binary systems exhibit a broad range of
variability that is not intrinsic to the individual components, we exclude
these stars from our classifier.

Figure 3 shows the makeup of our sample. Approximately 30% of our sample (2550
stars) belong to the Miscellaneous Variable and Unknown/Candidate classes,
which we do not use to train our classifier; instead, we use the classifier to
assign tentative classifications in §4. The rest of the sample is dominated by
luminous OBA stars and cool supergiants, with very few LBVs, OB[e] stars, and
WR. This is unsurprising given these stars’ high luminosity in the Gaia
bandpass, and the expected lifetimes of these evolutionary phases relative to
the lifetimes of exotic emission line objects (Ekström et al., 2012). The
imbalances in available training data across different classes, along with the
extreme sparsity of training data in the rare classes, will impact the
performance of the classifier if not properly addressed. (Chawla, 2010). We
discuss this issue for this particular sample in §4.1.

Figure 5: $M_{G}$ vs. $G-J$ for putative massive stars. Left: Stars are
colored by their label. Right: Stars are colored by their coarse label.
Contours for each coarse class correspond to 0.5 and 0.1 times the maximum
value of a kernel density estimate of the distribution of each class in the
CMD. Both panels illustrate a key difficulty faced by our classifier: both the
refined and coarse labels have significant overlap with each other in the CMD.
For example, the coolest/warmest YSGs have identical optical photometry to
RSGs/OBA supergiants respectively. Even in the coarse labels, the contours for
hot stars and emission line stars are nearly identical.

The left panel of Figure 5 shows the $M_{G}$ vs. $G-J$ CMD for all stars in
our sample that aren’t labeled as miscellaneous variables or
unknown/candidate, colored by their label. $G-J$ correlates reasonably well
with effective temperature in main sequence stars (Davenport & Covey, 2018).
From this plot, it is clear that many stars are misclassified in SIMBAD (the
worst example is one particular red star classified as an OBA star), reducing
the effectiveness of any machine learning algorithm, and propagating biases
into the results. Yellow supergiants are especially prone to this problem:
81/212 YSGs in the sample (38%) have $G-J<1$, consistent with the optical
colors of much hotter stars. Indeed, this problem would have been worse had we
not corrected for the presence of OBA stars in the initial sample of YSGs.
This issue may originate from bad distance estimates for individual stars
(which explains why our sample includes F and G dwarfs), bad estimates of
reddening (given our usage of monochromatic extinction coefficients),
previously unidentified variability, or the fact that many of these spectral
types were determined via stellar spectra taken on photographic plates (for
example, many spectral types for stars in the LMC come from Ardeberg et al.,
1972). We expect this issue to propagate in our results, increasing the
confusion between YSGs and hot stars.

Because we expect objects in some classes — especially those with an
evolutionary link such as main sequence, evolved, and supergiant OBA stars —
to appear similar in the training data set, we also assign all stars a coarse
label: all classes of OBA stars excluding OB[e] and OBAe are labeled “Hot”;
RSGs and YSGs are labeled “Cool”; WRs, LBVs, and both OB[e] and OBAe stars are
labeled “Emission” (EM for short); C/S/Giant stars and Yellow Dwarfs are
labeled “Contaminant”; and miscellaneous variables and unknown/candidates are
labeled “Unknown/Candidate.” The results of this labeling scheme are
summarized in Table 2 which shows the number of stars with a given refined
class that are assigned a particular coarse label. The right panel of Figure 5
shows the same CMD, with points colored by their coarse label. This leads to
some improvement: each coarse class lies in the approximate region of the CMD
that one would expect. Regardless, it is evident that selecting any one of
these classes solely from this optical photometry would be difficult: the
“cool” class has significant overlap with the “hot” class — largely driven by
the YSGs — emission line stars can be found at a range of $G_{BP}-G_{RP}$
colors, and there is significant overlap between low-mass contaminants that
will end their lives as white dwarves and true massive stars that will end
their lives in supernovae explosions. This point is emphasized by the
contours, which correspond to 0.5 and 0.1 times the maximum value of a kernel
density estimate of the distribution of each class.

We will use these coarse labels to train a second classifier. While these
coarse labels lose some specificity, each coarse class contains more stars,
hopefully increasing the performance of a classifier trained on these labels.
Furthermore, they still retain physical information while increasing the
number of stars in each class: the “cool” label contains stars with convective
envelopes, while “hot” stars contain radiative envelopes. Meanwhile, emission
line stars are notable for their variability. It is our hope that this second
classifier will still address two of our stated goals: to identify emission
line stars, and to reject contaminating low mass stars.

Table 2: Number of stars in a class that are assigned a given coarse label; we do not show with the Miscellanous Variable or Unknown/Candidate labels. | Coarse Label
---|---
Refined Label | Hot | Emission | Cool | Contaminant
Main Sequence OBA | 187 |  |  |
Evolved OBA | 409 |  |  |
Supergiant OBA | 798 |  |  |
OBA | 915 |  |  |
OBAe |  | 383 |  |
OB[e] |  | 12 |  |
WR |  | 37 |  |
LBV |  | 8 |  |
YSG |  |  | 212 |
RSG |  |  | 847 |
C/S/Giant |  |  |  | 118
Yellow Dwarf |  |  |  | 8
Total | 2309 | 440 | 1059 | 126

##  3 WISE Lightcurves

Variability in evolved massive stars has been well-characterized at timescales
from minutes to decades (e.g. Conroy et al., 2018; Dorn-Wallenstein et al.,
2019; Soraisam et al., 2020). In a study of massive stars in the Whirlpool
Galaxy (M51), Conroy et al. (2018) found that almost half of the stars
brighter than $M_{I}=-7$ were variable, with red stars nearing a variability
fraction of 1. Both red and extremely luminous blue stars exhibited quite high
amplitude ($\Delta I\geq 0.3$) variability. For spectral energy distributions
(SEDs) dominated by purely stellar light, mid-infrared (MIR) flux measurements
(and thus variability) is sensitive to (variations in) the bolometric
luminosity. However, for stars with significant circumstellar dust components
in their SEDs, MIR variability is correlated with both intrinsic bolometric
variability, and with dust creation/destruction processes in the circumstellar
medium (for example, in RSGs where it is correlated with the mass loss rate,
e.g. Yang et al., 2018).

The WISE mission provides lightcurves from stars in all parts of the sky,
observed over a $\sim$7 year baseline. Due to the scanning law adopted by WISE
(Wright et al., 2010), most stars not on the ecliptic poles are visited
approximately every $\sim 180$ days. All stars have a $\sim 3$-year data gap
from when WISE was placed in hibernation in February 2011 and when it was
reactivated in December 2013. WISE initially observed simultaneously in four
filters during its primary mission: $W1$ (3.4 $\mu$m), $W2$ (4.6 $\mu$m), $W3$
(12 $\mu$m), and $W4$ (22 $\mu$m). However, it was reduced to using only the
two bluest bands in its post-cryogenic survey mode called “NEOWISE”. The time,
duration, and number of individual observations during each $\sim 180$ day
visit depends on spatial geometry of the WISE scanning program, i.e. stars
closer to the ecliptic poles have longer duration visits (often exceeding a
week) with many epochs per visit, while star near the equator have very short
visits (typically a couple days) with only a few epochs per visit. Because
WISE lightcurves possess such non-uniform cadence, extracting detailed physics
for most individual stars is difficult. However, the WISE lightcurves place
fantastic constraints on MIR variability amplitude on longer timescales,
especially for evolved massive stars whose highest-amplitude variability
occurs over $\sim$year timescales. Such amplitude and timescale estimates are
related to the physical parameters of the star, potentially aiding in
classification.

For every star selected in §2 we queried the Single-Exposure (“L1b”) source
databases for all phases of the WISE mission, including the original 4-band,
partial cryogenic 3-band, and post-cryogenic 2-band NEOWISE tables. We used
astroquery to pull data in the region within 3 arcsec of the known source
location. To ensure high quality data for all recovered epochs, we require the
photometric quality flag to be PH_QUAL=A, the contamination flag to be
CC_FLAGS=00, the number of deblended sources flag to be NB=1, and the PSF
photometry fit quality (defined as the reduced $\chi^{2}$) in $W1$ to be
w1rchi2$<$5.

Only two of our stars did not have usable data from WISE: WISE
J074911.48-102000.2 (HD 63554) has no lightcurve available online, and WISE
J050128.62-701120.2, which does not have any corresponding object nearby on
SIMBAD. When calculating each of the variability metrics below, we instead
record a value of NaN (i.e., missing data). For the remaining stars, we ignore
the $W3$ and $W4$ data here due to the lower signal to noise and significantly
shorter observing baselines due to the loss of cryogenic observations after
the original WISE mission. As WISE observes simultaneously in all bands, we
can construct $W1-W2$ lightcurves without any interpolation, and simply
subtract the $W2$ data from the $W1$ data to obtain the $W1-W2$ color curve.
The left panels of Figure 6 show an example set of lightcurves for WISE
J000536.97+432405.0 (=HD 73), a B1.5IV star that illustrates the typical
observing cadence and variability of a bright star in our sample.

Figure 6: Example lightcurve for WISE J000536.97+432405.0. Top left: Raw
lightcurve, with $W1$ points plotted as blue errorbars, and $W2$ points
plotted in orange. Bottom left: Variability in $W1-W2$ plotted as green
errorbars. Top right: Binned $W1$ lightcurve. Blue points are binned data
(errorbars are smaller than the points). Black dotted line is the B-spline
interpolation. Time has been adjusted so the lightcurve is centered on $t=0$.
Bottom right: First-order derivative of the interpolant. Vertical blue lines
show the times where the derivative crosses zero, indicated by the horizontal
blue line.

###  3.1 Variability Metrics

####  3.1.1 Amplitude

For each of the three lightcurves of each object, we wish to extract simple
metrics that describe the amplitude and timescale of variability. We choose
$\chi^{2}$ about the median defined as

$\chi^{2}=\sum\big{(}\frac{M_{i}-\tilde{M}}{\sigma_{i}}\big{)}^{2}$ (2)

and the reduced-$\chi^{2}$

$\chi^{2}_{red}=\chi^{2}/(N-1)$ (3)

Where $M_{i}$ is a magnitude measurement, $\tilde{M}$ is the median of the
lightcurve, $\sigma_{i}$ is the corresponding error on the data point, and $N$
is the number of points in the lightcurve. We also calculate the Median
Absolute Deviation (MAD), and Error-Weighted MAD (EWM):

$\displaystyle MAD$ $\displaystyle=\mathrm{Median}(|M_{i}-\tilde{M}|)$ (4)
$\displaystyle EWM$
$\displaystyle=\mathrm{Median}(|M_{i}-\tilde{M}|/\sigma_{i})$ (5)

If the filtered and cleaned lightcurve only contains one good measurement (or
no good measurements), we automatically give it $\chi^{2}=\chi^{2}_{red}={\rm
MAD}={\rm EWM}={\tt NaN}$. We describe our method for treating missing data
below.

Figure 7: Top row: Distribution of derived $\chi^{2}_{red}$ (left) and $EWM$
(right) values in the $W1$ (blue), $W2$ (orange), and $W1-W2$ (green)
lightcurves. Bottom row: Scatter plots comparing different amplitude metrics.
The left panel shows $EWM$ vs. $\chi^{2}_{red}$ in the $W1$, $W2$, and $W1-W2$
lightcurves (using the same color-coding). The right panel shows
$\chi^{2}_{red}$ in $W2$ vs. $\chi^{2}_{red}$ in $W1$, with each point colored
by its coarse label.

The top panels of Figure 7 show the distributions of $\chi^{2}_{red}$ and EWM
derived for our sample. Values from $W1$ lightcurves are in blue, $W2$ in
orange, and $W1-W2$ in green. The bottom left panel shows a scatter plot of
$\chi^{2}_{red}$ vs. EWM for all three lightcurves. While the two measures
correlate reasonably well with each other, there is a branch of stars whose
lightcurves have high $\chi^{2}_{red}$ and low EWM; because the EWM is robust
to outliers, $\chi^{2}_{red}$ is an effective probe of lightcurves with sudden
brightening/fading events, while EWM is an effective selector for lightcurves
that display consistent variability. The bottom right panel shows a scatter
plot of $\chi^{2}_{red}$ in $W1$ vs. in $W2$, with each point colored by its
coarse label. A distinct branch of stars that are much more variable in $W2$
than $W1$ is clearly evident; oddly, the distributions of classes, EWM, and
broadband colors in this branch are consistent with the whole sample, and no
similar branch exists in the measured EWM values.

Visual inspection of the lightcurves of stars with $\chi^{2}_{red}<10$ in $W1$
and $\chi^{2}_{red}>100$ in $W2$ shows that these stars appear to have higher
signal to noise $W2$ measurements than $W1$, and have one observation during
which the star apparently becomes considerably redder, achieving $W1-W2$
values as high as $\sim 4$. Examining the times at which these extreme
reddening events occur shows a preference for times during the cryogenic WISE
survey, implying that this behavior is likely instrumental in origin, despite
our filtering using the provided quality flags. Nonetheless, we include
$\chi^{2}_{red}$ as it does not map perfectly onto $EWM$, and only 98 stars
fall in this regime. We do not yet know whether $\chi^{2}_{red}$, $EWM$, both
metrics, or neither are useful features for classification, so we keep both
with the intent of exploring their importance below.

####  3.1.2 Timescale

Many methods exist for estimating dominant timescales in lightcurves. Conroy
et al. (2018) use the Lomb-Scargle Periodogram (Lomb, 1976; Scargle, 1982) to
search for periodic variables. However, this approach suffers from numerous,
well-known issues (including accurate period recovery at low signal to noise),
and false peaks can easily be mistaken for real timescales, especially in
highly-irregularly sampled data, as is the case for the WISE lightcurves.
Soraisam et al. (2020) use a Gaussian process interpolation scheme coupled
with a wavelet analysis to estimate timescales in massive stars in M31
observed by the Palomar Transient Factory (PTF). However, with so few data
points, we found it difficult to obtain a reliable fit with a Gaussian
process, and even when the fit was successful, the resulting interpolant had a
large standard deviation in between WISE visits. The resulting measurements of
the characteristic timescale was more reflective of the kernel used.

Instead, we adopt a spline-based interpolation method, as follows. We first
subtract half of the sum of the times of the first and last available
observations, so that the lightcurve is centered at $t=0$. We then bin the
observations in each visit. Visits are defined as sets of points separated in
time by less than a defined threshold. Due to the WISE scanning law, some
stars near the ecliptic poles have visits separated by less than the typical
$\sim 180$ days. Therefore, we adopt 50 days as the threshold for visits. Two
stars in our sample are close enough to the ecliptic pole to be observed
nearly continuously such that we erroneously record two “visits”: one each
during the cryogenic and post-cryogenic surveys. However, neither star is
strongly variable and thus this small edge case does not substantially impact
our subsequent analyses.

For all observations in a given visit, we calculate the mean time and
$W1$/$W2$/$W1-W2$ measurement. We use scipy.interpolate.splrep in Python to
find the 3rd-order B-spline representation of the binned lightcurves, adopting
a smoothing factor s=10. This returns the knots, B-spline coefficients, and
degree of the spline. By definition, 3rd-order splines are differentiable, so
we use scipy.interpolate.splev to evaluate the first derivative of the spline
interpolant, and find the times when the derivative changes sign — i.e., when
the lightcurve reaches a maximum or minimum. As metrics of the characteristic
timescale of the lightcurve, we calculate the frequency of zero-crossings of
the first derivative of the spline interpolant, $\nu_{0}$ — calculated as the
number of times the derivative passes through zero, divided by the time
baseline of the lightcurve — $\langle\Delta t\rangle$, the mean of the
differences between successive zero-crossings, and the standard deviation of
the differences between successive zero-crossings, $\sigma_{\Delta t}$. For
stars with fewer than four visits, we automatically assign
$\nu_{0}=\langle\Delta t\rangle=\sigma_{\Delta t}={\tt NaN}$. The right panels
of Figure 6 show this process on the $W1$ lightcurve plotted in the left
panels. The blue points are the binned $W1$ measurements (the errors are
smaller than the size of the points), and the dotted black lines are the
spline interpolant (top right) and corresponding derivative (bottom right).

##  4 Machine Learning

###  4.1 Classifier Selection

The problem of classification based on broad-band photometry has a rich
history in the literature. With the advent of large surveys like the Sloan
Digital Sky Survey (SDSS, York et al. 2000), optical data could be coupled
with space-based MIR data to find the stellar locus in a 10-dimensional color-
space (Davenport et al., 2014). Recent efforts to separate stars from quasars,
or perform a regression on effective temperature with machine learning on
photometric data have been successful (Makhija et al., 2019; Bai et al.,
2019); however, these studies are often focused on main sequence, low mass
stars. This is an understandable choice given the rarity of evolved, high mass
stars, the absence of reliable distances to calculate luminosities from which
to select putative massive stars, and the fact that follow-up spectroscopy is
necessary in order to confirm a star’s membership in many important classes.

With the advent of Gaia DR2, luminosities can be easily determined, and
putative massive stars can be identified, as we do in §2. We wish to train an
algorithm that takes as input the broadband photometry and variability metrics
derived for our sample, and outputs spectral type classifications. Many
machine learning classifiers exist; of these, we wish to choose a flexible
model with well-understood mathematics, while avoiding techniques like neural
networks that can be difficult to interpret. Of the classifiers available in
the sklearn package, we decided to test a Random Forest (RF) classifier
(Breiman, 2001) — which consists of a collection of decision trees trained on
random subsets of samples and features — a Support Vector Machine (SVM)
classifier (Cortes & Vapnik, 1995) — which identifies hyperplanes in the
feature space that separate different classes — and a Gaussian process (GP)
classifier (Rasmussen & Williams, 2006). We refer the reader to these
publications, as well as to the sklearn documentation222https://scikit-
learn.org/stable/index.html for the mathematics and implementation details of
each classifier. In the multi-class case, a collection of classifiers are
trained on each possible pair of classes (one-versus-one or “ovo”), generating
a total of $N_{classes}(N_{classes}-1)/2$ classifiers where $N_{classes}$ is
the number of classes. Labels are assigned to test samples by allowing each
classifier to vote, and the label with the most votes is chosen (Knerr et al.,
1990).

Table 3: List of features passed to our machine learning classifiers, as well as clarifying definitions where relevant. WISE photometry used to calculate colors and magnitudes is from the ALLWISE data release (Cutri et al., 2013). All variability metrics are calculated from the WISE $W1$, $W2$, and $W1-W2$ lightcurves.  Feature | Definition
---|---
Colors & Magnitudes
$M_{G}$ | Absolute magnitude in Gaia $G$ band.
$G-J$ | From Gaia and 2MASS photometry.
$J-H$ | From 2MASS photometry.
$H-K_{s}$ | From 2MASS photometry.
$K_{s}-W1$ | From 2MASS and WISE photometry.
$W1-W2$ | From WISE photometry.
$W2-W3$ | From WISE photometry.
$W3-W4$ | From WISE photometry.
$M_{W1}$ | Absolute magnitude in WISE $W1$ band.
Variability Metrics
$\log\chi^{2}_{red}$ | Log of the reduced $\chi^{2}$.
$\log$ EWM | Log of the error-weighted Median Absolute Deviation.
$\nu_{0}$ | Frequency of zero-crossings of the first derivative of the
| spline interpolant.
$\log\langle\Delta t\rangle$ | Log of the average time between zero-crossings.
$\log\sigma_{\Delta t}$ | Log of the standard deviation of zero-crossing times.

Each type of classifier has a number of hyperparameters that affect the
performance of the classifier. For the RF classifier, n_estimators specifies
the number of trees in the forest, max_depth specifies how many branches each
decision tree in the forest can have, and max_features specifies the maximum
number of features each tree is trained on. We also set class_weight=balanced,
which weighs samples when fitting to account for the different frequencies of
each class in the data.

For the SVM classifier (SVC), $C$ is a regularization parameter that governs
the tradeoff between maximizing the margin and misclassifications in the
training set. Higher values of $C$ will force the SVC to correctly classify
every point, resulting in poor generalization (i.e. overfitting). The SVC
requires that the distance between two points in the feature space is defined
as the inner product of two vectors in the feature space,
$\langle\vec{X_{i}},\vec{X_{j}}\rangle$. Because the boundaries between
classes in our sample are not guaranteed to be linear, one can project the
samples into a much higher dimension space via a mapping function, $\Phi$,
where distances between two vectors in this space are calculated as
$\langle\vec{Z_{i}},\vec{Z_{j}}\rangle=\langle\Phi(\vec{X_{i}}),\Phi(\vec{X_{j}})\rangle$.
In reality, the transformed feature space can be incredibly high dimensional,
and explicitly mapping the data into this high-dimensional space is
computationally inefficient. Instead, we can adopt a kernel function, $K$ that
defines distances in the higher dimensional space, e.g.,
$K(\vec{X_{i}},\vec{X_{j}})=\langle\vec{Z_{i}},\vec{Z_{j}}\rangle$. Because
the kernel only takes the measured features as input and outputs a number,
using a kernel function implicitly maps the input feature space into a high
dimensional space without specifying $\Phi$.

Common choices of the kernel function include a linear kernel (i.e., the
Euclidean distance between individual samples in the feature space) and the
“radial basis function” (RBF) kernel:

$K(\vec{X_{i}},\vec{X_{j}})=e^{-\gamma||\vec{X_{i}}-\vec{X_{j}}||^{2}}$ (6)

where $\gamma$ governs the influence of the kernel function; lower values
result in increasingly linear boundaries, while high values result in the
decision function being entirely depended on individual points, creating small
islands of a given class centered on each training point. The advantage of
non-linear models like a SVM with a RBF kernel over linear methods is that the
decision boundaries can be much more flexible; the tradeoff is that the
contribution of individual features to the classifier cannot be easily
calculated without the unknown function $\Phi$ (see §4.3.1). The “optimum”
kernel and hyperparameters are chosen via a cross-validation strategy
described below. Finally, we also set class_weight=balanced for the SVC, which
automatically sets the value $C$ for class $i$ to $CN_{\rm samples}/(N_{\rm
classes}N_{i})$ where $N_{\rm samples}$ is the size of the sample, $N_{\rm
classes}$ is the number of classes, and $N_{i}$ is the number of objects in
the sample belonging to the class. This serves to weight rarer classes more
heavily. The GP classifier’s only hyperparameter is a choice of kernel, which
defines the covariance function of the Gaussian process.

To fit our classifiers, we first remove all stars labeled as miscellaneous
variables or unknown/candidates, as well as known binaries, and one star with
bad $J$ photometry, the Be star HD 53032. For features, we use the intrinsic
calculated value of $M_{G}$, as well as (uncorrected for extinction) $G-J$,
$J-H$, $H-K_{s}$, $K_{s}-W1$, $W1-W2$, $W3-W4$, and $\chi^{2}_{red}$, EWM,
$\nu_{0}$, $\langle\Delta t\rangle$, and $\sigma_{\Delta t}$ in all three
lightcurves — we indicate the WISE band or color that each variability metric
corresponds to with a subscript hereafter. The input features and a brief
description where relevant are listed in Table 3. Because $\chi^{2}_{red}$,
EWM, $\langle\Delta t\rangle$, and $\sigma_{\Delta t}$ have significant
dynamic range, we use the base-10 logarithm of these features. Features values
as well as labels for each star in our sample are given in Table 4.

Table 4: Feature values and assigned labels for all stars in our sample,
ordered by Right Ascension. Missing numbers are indicated with “-”.

Common Name |  $M_{G}$ [mag] |  $G-J$ [mag] |  $W1-W2$ [mag] | $\log\chi^{2}_{red,W1}$ |  $\log\langle\Delta t\rangle_{W1}$ [d] | Label | Coarse Label
---|---|---|---|---|---|---|---
HD 236270 | -3.54342624 | 0.3541 | 0.1420 | 0.389 | - | RSG | Cool
LS I +64 10 | -3.40866054 | 0.7154 | 0.0020 | - | - | OBA | Hot
LS I +60 69 | -3.22094985 | 0.8769 | -0.0340 | 0.061 | - | OBAe | EM
BD+62 2353 | -4.27358395 | 0.4653 | -0.0350 | 0.110 | 2.816 | SupergiantOBA | Hot
HD 73 | -3.29367442 | -0.6567 | -0.0520 | 3.157 | 2.701 | EvolvedOBA | Hot
HD 240496 | -3.84380534 | 1.0127 | 0.0130 | -0.064 | 2.764 | OBA | Hot
WISE J000559.28-790653.3 | -5.32554735 | 1.3391 | -0.0390 | 0.072 | 3.162 | Unknown/Candidate | Unknown/Candidate
LS I +59 30 | -3.34464936 | 0.6748 | -0.0430 | 0.234 | - | OBA | Hot
BD+57 2870 | -4.53103041 | 1.1711 | 0.0040 | 0.046 | 3.135 | SupergiantOBA | Hot
BD+62 1 | -3.31029252 | 0.9010 | 0.2660 | 1.622 | 2.813 | EvolvedOBA | Hot

Note. — This table is published in its entirety in the machine-readable
format. A portion is shown here for guidance regarding its form and content.

We now randomly split our sample into a training set with 70% of the samples,
and a test set with the remaining 30%, using a stratification strategy to
ensure the proportions of the classes in both sets are equal. The test set is
withheld until we are ready to assess the performance of the chosen
classifier. We then use sklearn.preprocessing.StandardScaler in Python to
scale the training data such that each feature has 0 mean and unit variance.
Because the data have missing values, we then use
sklearn.impute.IterativeImputer which uses a Bayesian ridge regression to
predict and replace missing values.

To test the accuracy of the imputer, we select only the rows from the training
set with no missing data. For each feature with missing data
($\log\chi^{2}_{red,W1-W2}$, $\log{\mathrm{E}WM}_{W1-W2}$, $\log\sigma_{\Delta
t,W1}$, $\log\sigma_{\Delta t,W2}$, and $\log\sigma_{\Delta t,W1-W2}$) we
randomly choose 200 objects, replace the value of the feature with NaN for
only these objects, transform the data using the scaler and imputer, and
calculate the fractional error between the true value and the imputed value.
The returned fractional errors for each feature centered around 0, and had a
low scatter with the exception of $\log{\mathrm{E}WM}_{W1-W2}$. However, this
has little impact on the classifier, as only two objects in the actual
training set have missing values for $\log{\mathrm{E}WM}_{W1-W2}$. We also
repeated this procedure for coarse labels: we select 200 random objects with a
given coarse label, replace a random feature from the list of features with
missing data with NaN for each object, and again calculate the fractional
error between the true and imputed values. We find that the imputer performs
poorly on Cool and Contaminant stars. Given that all of the features with
missing data are linked to variability, a significant fraction of red
supergiants display high amplitude variability (Conroy et al., 2018), the MIR
variability of AGB stars (which make up the bulk of the Contaminants) is
higher than RSGs in a given magnitude range (Yang et al., 2018), and our
sample of Cool and Contaminant stars contains objects in quite different
evolutionary states that nevertheless have similar colors and magnitudes, it
is unsurprising that the imputer is unable to predict the variability
properties of these stars. In §4.3.1, we will discuss the impact of these
features on the overall performance of our classifier.

For each classifier, we then initialize a corresponding sklearn classifier
object (e.g., sklearn.svm.SVC). To settle on the best values for the
hyperparameters we use sklearn.model_selection.GridSearchCV to perform a
cross-validation search on a grid of hyperparameters, using a stratified
K-fold strategy with $k=5$ to ensure that each fold has a representative
distribution of classes. For the RF, we search for n_estimators between 10 and
150 in steps of 10, max_depth between 10 and 100 in steps of 10, and allow
max_features to be either sqrt, log2, or None (where the maximum number of
features individual trees are trained on is the square root of, base-2
logarithm of, or equal to the number of features, see the documentation for
details). We also allow max_depth to take on the default value (None), such
that individual trees can be grown until each leaf only contains one sample.

For the SVC, we search for values of $C$ on a logarithmic grid with 1 dex
spacing between 0.01 and 100, and for the RBF kernel, we search for values of
$\gamma$ on a similar grid between 0.01 and 10. Additionally, we allow
$\gamma$ to be the default values of 1/n_features (where n_features is the
number of features). For the GP classifier, we only vary the kernel, as the
GaussianProcessClassifier object automatically optimizes the kernel
hyperparameters. We let the kernel be either linear, RBF, or the default (a
special case of the RBF kernel with the length scale equal to 1).

Each classifier object has a default method to score each set of
hyperparameters, e.g. the accuracy of predicted labels compared to true
labels. However, the classes in our training set are unbalanced (e.g. Figure
3), so inaccurately classifying every single LBV, for example, would have
little impact on the overall accuracy of the classifier. To account for this,
we instead use the balanced accuracy (Mosley, 2013; Guyon et al., 2015), which
weighs each sample by the frequency of that sample’s class in the training
set. Other options for scoring criteria exist, including some that help
maximize the classifier’s precision such as the weighted $F_{1}$ score and
Cohen’s kappa (Cohen, 1960). We experimented with using these scores, and
found that using the balanced accuracy minimizes misclassifications across all
classes (reflected in the diagonal in the left panels of Figures 9 and 12).
Note that this choice implicitly selects a classifier that performs well
across all classes, and is not optimized for specific classes. Future work
will explore the possibility of tuning a classifier to find specific classes
of rare stars.

Finally, we explore three variations of a voting classifier. Such a classifier
consists of an ensemble of individual classifiers, each which “votes” by
assigning a class to a given sample. The final assigned class can either be
chosen with a “hard” (the class with the most votes wins) or “soft” (class
assignments are weighted by the probabilities output by each classifier)
strategy. We construct two voting classifiers that each use a different voting
strategy, using RF, SVC, and GP classifiers as the individual components. We
refer to these as the Voting (Hard) and Voting (Soft) classifiers. We also
make a third voting classifier that also uses a soft voting strategy, but the
votes from each component classifier weighted by the balanced accuracy
determined via cross-validation. We refer to this as the Voting (Weighted)
classifier. We score each voting classifier by averaging the balanced accuracy
taken from five stratified folds of the data. Figure 8 shows the balanced
accuracy for the three optimized classifiers, as well as the three voting
classifiers; the SVC performs “best.” Both the Voting (Soft) and Voting
(Weighted) classifiers perform comparably with the worst classifier, the GP.
This is due to the fact that, while the SVC often selects one individual class
with high probability, both the RF and GP tend to select multiple classes with
high probability (with the GP sometimes selecting all classes with roughly
equal probability, usually slightly favoring the classes selected by the RF).
This can result in both the RF and GP voting for the wrong class with higher
probability than the correct vote from the SVC, leading to the poor observed
performance.

Figure 8: Balanced accuracy for each optimized classifier, averaged over five
foldings of the data. The three voting classifiers perform roughly
identically, with the Voting (Hard) classifier performing best.

###  4.2 SVC Performance

Figure 9: Left: Matrix showing the number of stars in the test set with true
label indicated on the y-axis that are assigned the label on the x-axis.
Center: Confusion matrix for the SVC, calculated by normalizing each row of
the left panel by the total number of stars in that row. Values correspond to
the fraction of samples in the test set with true label indicated on the
y-axis that are assigned the label on the x-axis, such that the values along
the diagonal are the fraction of each class that is correctly classified.
Right: Efficiency matrix for the SVC, calculated by normalizing each column of
the left panel by the total number of stars in that column. Values in each box
correspond to the fraction of samples in the test set assigned the label on
the x-axis that belong to the class on the y-axis, such that the values along
the diagonal correspond to the precision (one minus the contamination). Darker
colors in all panels correspond to more/a higher fraction of stars.

The procedure above results in values for the SVC hyperparameters of kernel =
linear, and C = 0.01. With these hyperparameters, we fit the SVC to the
training set, use the StandardScaler and IterativeImputer that were previously
fit to the training set to transform the test set, and use the SVC to predict
the labels of the test set. The left-hand panel of Figure 9 shows the raw
number of stars in the test set with the true label given on the x-axis, and
the predicted label given on the y-axis. The center and right-hand panels show
this matrix, where each row/column is normalized by the total number of stars
in that row/column, yielding the confusion/efficiency matrices, respectively.
The $i,j$ entry in the confusion matrix corresponds to the fraction of objects
in the test set belonging to class $i$ (shown on the y-axis) that are assigned
class $j$ (shown on the x-axis). Entries along the diagonal are the
completeness (also called the recall in some contexts), i.e., the percentage
of a given class that is accurately recovered by the classifier. The $i,j$
entry in the efficiency matrix is the fraction of objects in the test set
classified as $j$, that belong to class $i$. Entries along the diagonal are
equivalent to the precision (equivalent to one minus the contamination), i.e.,
the percentage of an observed class that is made up of true members of that
class. Figure 10 shows the completeness versus the contamination for each
class. Completeness is just the diagonal of the corresponding row/column of
the confusion matrix, and contamination is one minus the diagonal of the
corresponding row in the efficiency matrix.

Figure 10: Completeness versus contamination of each class in the test set. A
high completeness value implies members of that class are accurately
classified, while a low contamination value implies that an object classified
as such is likely to belong to that class. The figure is roughly divided into
four quadrants; stars with classes in the bottom-right quadrant can be
considered to be well-classified, in the sense that they have high
completeness and low contamination.

The SVC performs poorly on non-supergiant OBA stars. This is perhaps
unsurprising given that both the observed colors and interior structures of
OBA stars as they evolve from the zero-age main sequence (ZAMS) to the
terminal-age main sequence (TAMS) do not change drastically compared to the
much more evolved states that we also consider. The classifier classifies OBAe
stars with somewhat lower contamination compared to main sequence and evolved
OBA stars, though with comparably low completeness. True OBAe stars are either
misclassified as other types of OBA star or as WRs, while stars falsely
labeled as OBAe are mostly true OBA stars, with the exception of one true WR.
75% of WR stars are recovered, but only 6/30 stars identified as WRs in the
test set are true WRs; given the importance of WRs for both the physics of
mass loss and studying evolved massive stellar populations (Dorn-Wallenstein &
Levesque, 2018, 2020), future work will focus on developing a classifier
specifically for identifying WR stars.

All LBVs in the test set are recovered; while such high accuracy is often seen
as a sign of overfitting, we choose not to focus on this subclass, given both
the disputed evolutionary status of LBVs (Smith & Tombleson, 2015; Humphreys
et al., 2016; Aadland et al., 2018) and the fact that only two LBVs exist in
the training set. Yellow supergiants are only classified with 27% accuracy. As
discussed in §2 and shown in Figure 5, the yellow supergiant label is assigned
to stars with optical colors consistent with hot stars as well as RSGs. This
is reflected in the types of stars that YSGs are mistaken for, as well as the
stars that are mistaken for YSGs. Overall, the classifier performs best on the
coolest stars in the sample. RSGs are classified with 96% accuracy, and only
10% contamination. Meanwhile the classifier performs exceptionally well at
identifying low mass contaminants, at the cost of misclassifying four RSGs,
two OBAe stars, and one OBA star.

Overall, an SVC trained on these refined labels appears to have little use.
With the exception of RSGs and low-mass giants, the remaining classes have low
accuracy, high contamination, or both. We nonetheless use the SVC to predict
labels for the 2550 stars initially labeled as “Miscellaneous Variable” or
“Unknown/Candidate.” We identify 79 candidate RSGs and 36 candidate C/S/Giant
stars, of which we expect $\sim$71 and 30 to be genuine, respectively, given
the efficiency matrix. We list the candidate RSGs in Table 5. A small
spectroscopic observing campaign would be easily confirm the ability of this
classifier to correctly identify RSGs and low-mass giants.

Table 5: Common names and coordinates of stars predicted to be RSGs by the SVC
trained on refined labels.

Common Name | R.A. [deg] | Dec [deg]
---|---|---
SP77 48-11 | 81.07900941 | -70.43417562
WISE J185608.58-163255.1 | 284.03575762 | -16.54867009
W61 19-14 | 83.07777354 | -67.52941938
OGLE BRIGHT-LMC-MISC-169 | 72.94711333 | -69.32348227
WISE J194127.64+385155.3 | 295.36520609 | 38.86536427
NGC 2004 BBBC 431 | 82.69161818 | -67.29036242
[KWV2015] J045626.51-692350.6 | 74.11062177 | -69.39740804
W61 6-54 | 85.54018822 | -69.21978048
WISE J064232.30-715243.3 | 100.63462144 | -71.87871874
W61 6-34 | 85.51621031 | -69.21870016

Note. — This table is published in its entirety in the machine-readable
format. A portion is shown here for guidance regarding its form and content.

###  4.3 Performance on Coarse Labelling

Examining Figures 9 and 10, we see that, while the classifier is not
especially accurate except for the coolest stars, the classifier is roughly
useful for sorting the test set into broad categories: different types of OBA
stars are mostly (mis)classified as other classes of OBA stars; the same is
true for emission line stars (OBAe, OB[e], WR, and LBV) and cool stars (YSG,
RSG, C/S/Giant).

For this reason, we also utilize the coarse labels introduced in §2. We repeat
the entire process described above, beginning with the selection of the
classifier. Figure 11 shows the balanced accuracy for each of the classifiers
discussed above, trained on the coarse labels, using a five-fold cross-
validation to optimize the hyperparameters of each classifier. We find that
once again, the SVC yields the highest balanced accuracy (0.876).

Figure 11: Similar to Figure 8 for classifiers trained on the coarse labels.

We keep the same scaled and imputed training and testing sets, and perform a
five-fold cross-validation as before to find the optimal hyperparameters,
which are kernel=rbf, C = 1, and gamma = 1/n_features. With these
hyperparameters, we fit the SVC to the training set before predicting labels
for the test set. Figure 12 shows the confusion and efficiency matrices
similar to Figure 9, while Figure 13 shows the completeness versus the
contamination similar to Figure 10. All told, the SVC performs significantly
better compared to the classifier trained on the refined labels, recovering
all classes with $\geq 75$% completeness and $\leq 30$% contamination.

Of the emission line stars that are correctly identified, 83 are OBAe stars,
three are OB[e] stars, eight are WRs, and two are LBVs. This is 73%, 75%,
100%, and 100%, respectively of these stars that are in the test set, implying
that the performance of the SVC on emission line stars is not dominated
entirely by OBAe stars (which comprise the majority of emission line stars in
the test set). Of the stars mislabeled as contaminants, two are OBA stars and
two are RSGs. One true C/S/Giant star, and two Yellow Dwarfs are
misclassified.

Figure 12: Similar to Figure 9 for the SVC trained using the coarse labels.
Note that significantly more stars fall along the diagonal of each plot,
reflecting the improved performance of the SVC on the coarse labels. Figure
13: Similar to Figure 10 for the SVC trained using the coarse labels. All
coarse classes have high completeness and low contamination. Figure 14:
Distribution of coarse labels assigned to 2550 stars with no previously known
class.

We then use the SVC to predict the coarse labels for the same 2550 stars as
above. Figure 14 shows the distribution of these predicted labels. The
majority (2472 stars) are labeled as Hot. 63 stars are labeled as Cool, three
of which are already identified in SIMBAD as candidate AGBs or RGBs. 14 of
these stars are labeled as emission line stars, of which 9-10 are likely to
actually be emission line stars, assuming 30% contamination. We list all 2550
stars’ common names, coordinates, and predicted coarse label in Table 6.

####  4.3.1 Feature Importance

We can also identify which features contribute most to the overall performance
of the classifier on the coarse labels. To do this, we initialize a new SVC
object with the same hyperparameters, and perform a greedy search over
features. For each feature in the scaled and imputed training set, we train
the SVC on just this feature across five stratified folds of the training set,
and record the average and standard deviation of the balanced accuracy. We
select the feature that yields the highest average balanced accuracy. This has
the advantage of ensuring the contribution of each feature to the balanced
accuracy is stable across subsets of the data. We then train the SVC on all
combinations of this feature and the remaining features, selecting which
combination again yields the highest average balanced accuracy. This process
is repeated until all features are used.

Figure 15 illustrates this process. The x-axis shows the feature that is
selected at each stage of the greedy search. The y-axis shows the mean
balanced accuracy of the SVC at that stage. Errorbars show the standard
deviation of the balanced accuracy across the five folds of the training set.
The balanced accuracy reaches a maximum after the first seven features: $J-H$,
$W1-W2$, $M_{W1}$, $K_{s}-W1$, $\log EWM_{W2}$, $\log\chi^{2}_{red,W2}$, and
$\log EWM_{W1-W2}$. However, the contribution to the balanced accuracy from
all but the first four features is small. This suggests that variability
amplitude is a useful, though not critical metric to obtain, while variability
timescales are not necessary. Finally this suggests that photometry bluer than
$J$-band is also unnecessary.

Figure 15: Mean balanced accuracy of the SVC for coarse labels trained on
successively added features, calculated from five stratified folds of the
data. The balanced accuracy reaches a maximum after the first seven features.
Errorbars indicate the standard deviation of the balanced accuracy across
folds.

We can also examine the importance of each feature for classifying individual
classes. We perform the same greedy search over the features, instead
calculating a performance metric that focuses on the performance on a specific
class. One option is the $F_{\beta}$ measure:

$F_{\beta}=(1+\beta^{2})\frac{{\rm completeness}\cdot{\rm
precision}}{\beta^{2}\cdot{\rm completeness}+{\rm precision}}$ (7)

where $\beta$ is a free parameter that sets the relative importance of
completeness compared to precision. Common choices are $\beta=1$ (i.e.,
$F_{1}$, a harmonic mean of completeness and precision), $\beta=0.5$, and
$\beta=2$ (Chinchor & D, 1992). We adopt $F_{2}$ (i.e., $\beta=2$), because we
prioritize generating complete samples of rare massive stars.

The left panels of Figure 16 show the $F_{2}$ measure as a function of
successively added features, calculated specifically for hot stars (top),
emission line stars (second panel), cool stars (third panel), and contaminants
(bottom). The results for both hot and cool stars are mostly similar to the
results for the overall classifier in Figure 15, in the sense that the best
performance is reached after including a mix of near- and mid-infrared colors
and magnitudes. The main difference is that $M_{G}$ is the fourth most
important feature for classifying hot stars.

For emission line stars, a maximum in the mean $F_{2}$ is reached after 11
features: $W1-W2$, $M_{W1}$, $G_{J}$, $\log\langle\Delta t\rangle_{W1}$,
$\log\langle\Delta t\rangle_{W2}$, $\log\sigma_{\Delta t,W2}$, $W3-W4$,
$W2-W3$, $K_{s}-W1$, $\nu_{0,W2}$, and $J-H$. However, given the errorbars,
only the first three features contribute meaningfully, with the remaining
features consistent with a constant value of $F_{2}$. Interestingly, compared
to its contribution to the overall balanced accuracy of the classifier, bluer
photometry (signified by the presence of $G-J$ in the above list) is much more
important for identifying emission line stars. While variability metrics are
included in the above list, they do not significantly contribute to the
$F_{2}$ score.

For contaminants, a total of 15 features are required in order to maximize
$F_{2}$: $J-H$, $\log\sigma_{\Delta t,W1-W2}$, $\log\sigma_{\Delta t,W2}$,
$\log\langle\Delta t\rangle_{W1}$, $\log\langle\Delta t\rangle_{W1-W2}$,
$\nu_{0,W2}$, $\log\sigma_{\Delta t,W1}$, $\nu_{0,W1-W2}$, $M_{W1}$,
$H-K_{s}$, $K_{s}-W1$, $\log$ EWMW2, $\log\chi^{2}_{red,W1-W1}$, $\log$
EWMW1-W2, and $\log\chi^{2}_{red,W1}$. Notably, the $F_{2}$ measure first
decreases as features are added, before increasing to the maximum after
$M_{W1}$. This trend is unintuitive compared to the other panels in the
Figure. It may be a result of the fact that increased features improve the
precision of the classifier at the cost of completeness, resulting in a
decrease in $F_{2}$ due to the increased weighting of completeness.

Figure 16: Left panels are similar to Figure 15, except using the $F_{2}$
measure calculated for hot stars (top), emission line stars (second panel),
cool stars (third panel), and contaminants (bottom). The right panel in each
row shows a scatter plot of only first and second most important features
(indicated with blue and red text respectively) drawn from the test set, with
stars belonging to the corresponding class in each row highlighted. Note that
the features plotted are the scaled and imputed values, not the original
values listed in Table 4.

To demonstrate the capabilities of the classifier using a limited set of
features, we plot the scaled and imputed test set — which was not used in the
greedy search algorithm — in the right panels of Figure 16, using only the two
most important features in each row. Stars belonging to the corresponding
coarse class are plotted as larger, colored points. In all cases, most members
of the test set with that label are well-separated from the other stars. As
expected from the left panels, most of the separation is along the x-axis
which corresponds to the most important feature, with the second-most
important feature plotted on the y-axis providing some additional
differentiation, especially for hot stars. In the case of contaminants, the
second-most important feature provides little-to-no additional information,
consistent with the lack of change in $F_{2}$ with increased features.

We conclude that while small numbers of features can be used to identify hot,
cool, and (remarkably) emission line stars, a large number of features is
necessary in order to maximize the number of accurately identified
contaminants. This includes time domain features, where the drastically
different structures of old AGB and RGB stars compared to massive cool
supergiants may be imprinted. Already extragalactic massive star samples are
contaminated by foreground giants in the Milky Way halo; distant stars that
can be resolved by JWST and Roman will have comparable brightnesses to cool
dwarfs that are too faint to be filtered out using astrometry from Gaia .
Developing the infrastructure to reliably remove these contaminating objects
from massive star samples will be that much more critical.

##  5 Discussion & Conclusion

In the coming decades, space-based infrared observatories like Webb and Roman
will give us access to unprecedentedly large samples of evolved massive stars.
Therefore, we need to be prepared to leverage these data to search for stars
in the most interesting evolutionary states. Obtaining spectroscopy of
individual stars does not scale well at the size of the expected samples,
while linear cuts in color-magnitude space are too simplistic and ignore
emission line objects. Here we have demonstrated the promising performance of
a support vector machine trained on $\sim 0.5-22$ $\mu$m photometry and simple
variability metrics. However, with currently available labels, we are not able
to construct a classifier that performs well at the level of granularity
needed for many science cases.

Our main results are summarized as follows:

  * •

We have assembled a large sample of evolved massive stars using distances from
Gaia DR2 and Bailer-Jones et al. (2018), with high precision infrared
photometry from Gaia , 2MASS, and WISE.

  * •

Using SIMBAD, we assign labels to all stars, and find that the sample contains
a number of low mass contaminants.

  * •

We find that, of the classification methods we applied, a support vector
machine classification (SVC) algorithm is best at accurately labeling evolved
massive stars. The SVC is fast, and has the added benefit that the underlying
mathematics are well-understood.

  * •

The SVC trained on refined labels is capable of identifying low mass red giant
contaminants with high accuracy. However, the overall performance of this
classifier is quite poor, and we do not recommend its use at present.

  * •

The SVC trained using coarse spectral types performs better, as measured with
the balanced accuracy score. We find higher completeness and lower
contamination (Figures 12 and 13) compared to the SVC trained on the refined
labels (Figures 9 and 13). With this classifier, we identify 14 candidate
emission line stars from a sample of $\sim 2500$ unlabeled stars. We plan to
obtain spectroscopy of these stars to confirm our results.

  * •

We find that the SVC performs equally as well with only a small subset of
features. These features are mostly infrared colors and absolute magnitudes —
i.e., those least affected by reddening — with small contributions from
infrared variability metrics. However, if we change our performance metric to
one that focuses on emission line stars, optimal performance of the classifier
requires some red-optical photometry. We find that the added benefit of using
variability metrics may not be worth the investment in telescope time in order
to measure them.

Ultimately, the performance of the SVC trained on the refined labels is poor.
All stars in the sample are bright ($W1<14$), and the input features we use
are easily measured, implying we are not limited by the quality of the data.
However, the labelling itself is not of sufficient accuracy, as can be seen in
Figure 5. Labels are derived inhomogeneously, and many are from spectroscopy
that is now more than 50 years old. Unfortunately, these are the best labels
available for this sample. At the present, though curated lists of different
subclasses of massive stars exist (e.g. Richardson & Mehner, 2018), no unified
catalog of massive stars in our Galaxy or the Magellanic Clouds exists.

Modern all-sky surveys have already given us access to precision photometric
and spectroscopic measurements of unprecedented numbers of stars. Massive
stars are bright, and so the existing data is of suitable signal to noise to
perform spectroscopic classification. However, they are often excluded from
analyses that provide value-added measurements like effective temperatures,
surface gravities, compositions, radial velocities, and more that can be used
to accurately classify massive stars. In order to prepare ourselves for the
era of Webb and Roman, we must develop pipelines specifically tuned for
evolved massive stars. This is especially true for the classes that are
underrepresented in our dataset, i.e., rare emission line stars.

Along with better labels, more data will become available via future data
releases of the Gaia mission. The recent early third Gaia release contains
modest improvements in precision and sample size that are unlikely to affect
our results given the high quality of the photometry in our sample. However,
the full Gaia DR3 will contain low-resolution spectra as well as epoch
photometry for a limited number of sources, which have the potential to
significantly improve the performance of a machine learning classifier. On the
horizon, the Legacy Survey of Space and Time (LSST) conducted at the Vera
Rubin Observatory will measure the multi-color variability of massive stars at
higher cadence, while its large telescope aperture will help define a much
larger sample. As we demonstrate with Figure 16, it is possible to select
features that maximize the performance of the SVC for specific classes. With a
larger sample, we may be able to optimize the SVC to search for specific
classes of evolved massive stars.

The authors acknowledge that the work presented was largely conducted on the
traditional land of the first people of Seattle, the Duwamish People past and
present and honor with gratitude the land itself and the Duwamish Tribe.
This research was supported by NSF grant AST 1714285 awarded to EML. JRAD and
DH acknowledge support from the DiRAC Institute in the Department of Astronomy
at the University of Washington. The DiRAC Institute is supported through
generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences,
and the Washington Research Foundation. DH is supported by the Women In
Science Excel (WISE) programme of the Netherlands Organisation for Scientific
Research (NWO). This project was developed in part at the 2018 Gaia Sprint,
hosted by the eScience and DiRAC Institutes at the University of Washington,
Seattle. This research has made use of the VizieR catalogue access tool, CDS,
Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the
VizieR service was published in A&AS 143, 23. This research has made use of
the SIMBAD database, operated at CDS, Strasbourg, France. This publication
makes use of data products from the Two Micron All Sky Survey, which is a
joint project of the University of Massachusetts and the Infrared Processing
and Analysis Center/California Institute of Technology, funded by the National
Aeronautics and Space Administration and the National Science Foundation. This
work made use of the following software:

##  Appendix A Coarse labels for 2550 stars

Table 6 lists all 2550 stars with no known label, and predicted labels
generated by the SVC trained on coarse labels.

Table 6: Common names, coordinates, and predicted labels of 2550 stars input
to the SVC trained on coarse labels, ordered by Right Ascension.

Common Name | R.A. [deg] | Dec [deg] | Predicted Coarse Label
---|---|---|---
WISE J000559.28-790653.3 | 1.49713706 | -79.11483482 | Hot
TYC 4500-1480-1 | 2.86210879 | 79.08686958 | Hot
BD+61 45 | 5.25504760 | 62.77064970 | Hot
NGC 104 LEE 2520 | 5.41170226 | -72.21106679 | Hot
WISE J002203.44-693554.7 | 5.51434821 | -69.59851087 | Hot
WISE J002207.43-742212.1 | 5.53102165 | -74.37003199 | Hot
WISE J002318.05-742326.4 | 5.82523611 | -74.39068759 | Hot
WISE J002340.20-750446.9 | 5.91756693 | -75.07972556 | Hot
WISE J002758.92-764527.2 | 6.99552600 | -76.75757402 | Hot
WISE J002759.32-742119.8 | 6.99734043 | -74.35552728 | Hot

Note. — This table is published in its entirety in the machine-readable
format. A portion is shown here for guidance regarding its form and content.

## References

  * Aadland et al. (2018) Aadland, E., Massey, P., Neugent, K. F., & Drout, M. R. 2018, AJ, 156, 294, doi: 10.3847/1538-3881/aaeb96
  * Ardeberg et al. (1972) Ardeberg, A., Brunet, J. P., Maurice, E., & Prevot, L. 1972, Astronomy and Astrophysics Supplement Series, 6, 249
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bai et al. (2019) Bai, Y., Liu, J., Wang, S., & Yang, F. 2019, AJ, 157, 9, doi: 10.3847/1538-3881/aaf009
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58, doi: 10.3847/1538-3881/aacb21
  * Breiman (2001) Breiman, L. 2001, Mach. Learn., 45, 5–32, doi: 10.1023/A:1010933404324
  * Chawla (2010) Chawla, N. V. 2010, Data Mining for Imbalanced Datasets: An Overview, ed. O. Maimon & L. Rokach (Boston, MA: Springer US), 875–886.  https://doi.org/10.1007/978-0-387-09823-4_45
  * Chinchor & D (1992) Chinchor, N., & D, P. 1992, in In Proceedings of the Fourth Message Understanding Conference, 22–29
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Cohen (1960) Cohen, J. 1960, Educational and Psychological Measurement, 20, 37, doi: 10.1177/001316446002000104
  * Conroy et al. (2018) Conroy, C., Strader, J., van Dokkum, P., et al. 2018, ApJ, 864, 111, doi: 10.3847/1538-4357/aad460
  * Conti et al. (1983) Conti, P. S., Garmany, C. D., De Loore, C., & Vanbeveren, D. 1983, ApJ, 274, 302, doi: 10.1086/161447
  * Cortes & Vapnik (1995) Cortes, C., & Vapnik, V. 1995, Machine Learning, 20, 273, doi: 10.1007/BF00994018
  * Cutri et al. (2013) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, Explanatory Supplement to the AllWISE Data Release Products, Explanatory Supplement to the AllWISE Data Release Products
  * Davenport & Covey (2018) Davenport, J. R. A., & Covey, K. R. 2018, ApJ, 868, 151
  * Davenport et al. (2014) Davenport, J. R. A., Ivezić, Ž., Becker, A. C., et al. 2014, MNRAS, 440, 3430, doi: 10.1093/mnras/stu466
  * Dorn-Wallenstein & Levesque (2018) Dorn-Wallenstein, T. Z., & Levesque, E. M. 2018, ArXiv e-prints.  https://arxiv.org/abs/1810.01902
  * Dorn-Wallenstein & Levesque (2020) —. 2020, arXiv e-prints, arXiv:2004.13040.  https://arxiv.org/abs/2004.13040
  * Dorn-Wallenstein et al. (2019) Dorn-Wallenstein, T. Z., Levesque, E. M., & Davenport, J. R. A. 2019, ApJ, 878, 155, doi: 10.3847/1538-4357/ab223f
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146, doi: 10.1051/0004-6361/201117751
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Helmi, A., van Leeuwen, F., et al. 2018b, A&A, 616, A12, doi: 10.1051/0004-6361/201832698
  * García Pérez et al. (2016) García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2016, AJ, 151, 144, doi: 10.3847/0004-6256/151/6/144
  * Ginsburg et al. (2019) Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98, doi: 10.3847/1538-3881/aafc33
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279, doi: 10.1086/376774
  * Götberg et al. (2018) Götberg, Y., de Mink, S. E., Groh, J. H., et al. 2018, A&A, 615, A78, doi: 10.1051/0004-6361/201732274
  * Guyon et al. (2015) Guyon, I., Bennett, K., Cawley, G., et al. 2015, in 2015 International Joint Conference on Neural Networks (IJCNN), 1–8
  * Humphreys et al. (2016) Humphreys, R. M., Weis, K., Davidson, K., & Gordon, M. S. 2016, ApJ, 825, 64, doi: 10.3847/0004-637X/825/1/64
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Jones et al. (2017) Jones, O. C., Meixner, M., Justtanont, K., & Glasse, A. 2017, ApJ, 841, 15, doi: 10.3847/1538-4357/aa6bf6
  * Knerr et al. (1990) Knerr, S., Personnaz, L., & Dreyfus, G. 1990, Neurocomputing: Algorithms, Architectures and applications, 1, doi: 10.1007/978-3-642-76153-9_5
  * Kovács (2000a) Kovács, G. 2000a, A&A, 360, L1.  https://arxiv.org/abs/astro-ph/0007271
  * Kovács (2000b) —. 2000b, A&A, 363, L1.  https://arxiv.org/abs/astro-ph/0011056
  * Levesque et al. (2012) Levesque, E. M., Leitherer, C., Ekstrom, S., Meynet, G., & Schaerer, D. 2012, ApJ, 751, 67, doi: 10.1088/0004-637X/751/1/67
  * Levesque et al. (2005) Levesque, E. M., Massey, P., Olsen, K. A. G., et al. 2005, ApJ, 628, 973, doi: 10.1086/430901
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2, doi: 10.1051/0004-6361/201832727
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * Mainzer et al. (2011) Mainzer, A., Bauer, J., Grav, T., et al. 2011, ApJ, 731, 53, doi: 10.1088/0004-637X/731/1/53
  * Makhija et al. (2019) Makhija, S., Saha, S., Basak, S., & Das, M. 2019, Astronomy and Computing, 29, 100313, doi: 10.1016/j.ascom.2019.100313
  * Malhan et al. (2018) Malhan, K., Ibata, R. A., & Martin, N. F. 2018, MNRAS, 481, 3442, doi: 10.1093/mnras/sty2474
  * Massey et al. (2007) Massey, P., Olsen, K. A. G., Hodge, P. W., et al. 2007, AJ, 133, 2393, doi: 10.1086/513319
  * Massey et al. (2006) —. 2006, AJ, 131, 2478, doi: 10.1086/503256
  * Massey et al. (2009) Massey, P., Silva, D. R., Levesque, E. M., et al. 2009, ApJ, 703, 420, doi: 10.1088/0004-637X/703/1/420
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Mosley (2013) Mosley, L. 2013, in A balanced approach to the multi-class imbalance problem
  * Neugent et al. (2018a) Neugent, K. F., Levesque, E. M., & Massey, P. 2018a, AJ, 156, 225, doi: 10.3847/1538-3881/aae4e0
  * Neugent et al. (2020) Neugent, K. F., Massey, P., Georgy, C., et al. 2020, ApJ, 889, 44, doi: 10.3847/1538-4357/ab5ba0
  * Neugent et al. (2018b) Neugent, K. F., Massey, P., & Morrell, N. 2018b, ArXiv e-prints.  https://arxiv.org/abs/1807.01209
  * Ochsenbein et al. (2000) Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23, doi: 10.1051/aas:2000169
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Price-Whelan et al. (2018) Price-Whelan, A., Mechev, A., & jumeroag. 2018, adrn/makecite: v0.1, doi: 10.5281/zenodo.1343295.  https://doi.org/10.5281/zenodo.1343295
  * Rasmussen & Williams (2006) Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning
  * Richardson & Mehner (2018) Richardson, N. D., & Mehner, A. 2018, Research Notes of the American Astronomical Society, 2, 121, doi: 10.3847/2515-5172/aad1f3
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Smith & Tombleson (2015) Smith, N., & Tombleson, R. 2015, MNRAS, 447, 598, doi: 10.1093/mnras/stu2430
  * Soraisam et al. (2020) Soraisam, M. D., Bildsten, L., Drout, M. R., et al. 2020, ApJ, 893, 11, doi: 10.3847/1538-4357/ab7b7b
  * Spergel et al. (2013) Spergel, D., Gehrels, N., Breckinridge, J., et al. 2013, arXiv e-prints, arXiv:1305.5425.  https://arxiv.org/abs/1305.5425
  * Stanway et al. (2020) Stanway, E. R., Eldridge, J. J., & Chrimes, A. A. 2020, MNRAS, 497, 2201, doi: 10.1093/mnras/staa2089
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, ArXiv e-prints.  https://arxiv.org/abs/1801.02634
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, ArXiv e-prints, arXiv:1102.1523.  https://arxiv.org/abs/1102.1523
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9, doi: 10.1051/aas:2000332
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yang et al. (2018) Yang, M., Bonanos, A. Z., Jiang, B.-W., et al. 2018, A&A, 616, A175, doi: 10.1051/0004-6361/201832833
  * Yang et al. (2019) —. 2019, A&A, 629, A91, doi: 10.1051/0004-6361/201935916
  * York et al. (2000) York, D. G., Adelman, J., Anderson, Jr., J. E., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513