THE MILKY WAY PROJECT: LEVERAGING CITIZEN SCIENCE AND MACHINE LEARNING TO DETECT INTERSTELLAR BUBBLES

Generally speaking, a bubble is a shell-like, 1–30 parsec-scale cavity in the ISM, cleared by a combination of thermal overpressure, radiation pressure, and stellar winds. The basic structure of a bubble is shown in Figure 1. Strömgren (1939) first derived the ionization structure around OBA stars in the ISM. Stellar radiation both ionizes and heats gas up to ∼10⁴ K, creating a strong overpressure that drives expansion into the surrounding medium. This expansion creates an overdense shell of gas along the expansion front. Castor et al. (1975) and Weaver et al. (1977) extended this analysis to include the effects of strong stellar winds from O and early B stars. The main effect of a stellar wind is to create a shocked, high-temperature (10⁶ K) region within the 10⁴ K ionization region. This shock sweeps up additional material within the ionization region, potentially creating a second shell. Though the term "bubble" originally referred to cavities cleared by stellar winds, modern observational studies of "bubbles" tend to include both these objects and classical H ii regions. Furthermore, Beaumont & Williams (2010) demonstrated that many bubbles are embedded in relatively thin clouds, and more closely resemble rings than spheres. Throughout this paper, we use the term "bubble" broadly to refer to any ring- or shell-like cavity cleared by a young or main sequence star.

Zoom In Zoom Out Reset image size

The Spitzer Space Telescope and its surveys of the Galactic midplane—glimpse (Benjamin et al. 2003) and mipsgal (Carey et al. 2009)—enabled comprehensive, statistical studies of bubbles in the Galaxy. Mid-infrared wavelengths are well-suited for bubble observations, as they penetrate deeper into the Galactic disk and match bubble emission features. Figure 1 schematically depicts the observational signature of a bubble in Spitzer data. The interface between bubble shells and the ambient ISM excites polycyclic aromatic hydrocarbon (PAH) emission features, several of which fall within Spitzer's 8 μm bandpass. Bubble interiors often emit at 24 μm, due to emission from hot dust grains (Everett & Churchwell 2010).

Churchwell et al. (2006, 2007) carried out the first search for bubbles in Spitzer images of the Galactic plane, yielding a catalog of some 600 bubbles located at |ℓ| < 65°, |b| < 1°. These objects were identified by four astronomers manually searching through 3.5–8.0 μm images (they did not have access to 24 μm images). Churchwell et al. (2006) noted that each astronomer possessed different selection biases, and cautioned that their catalog was likely incomplete.

In an attempt to overcome the inherent bias in manual classification, the web-based Milky Way Project (Simpson et al. 2012) enlisted over 35,000 citizen scientists to search for bubbles in Spitzer data. The Milky Way Project (hereafter MWP) presented color-composite Spitzer images at 4.5 μm, 8 μm, and 24 μm, and asked citizen scientists to draw ellipses around potential bubbles. Figure 2 shows two typical images presented to the public. Each image in the MWP was 800 × 400 pixels, with a pixel scale ranging between 135 and 675 pixel⁻¹. The citizen science effort produced a dramatically larger catalog of ∼5000 objects, nearly 10 times the number in the catalog of ∼600 shells cataloged by the four astronomers of the Churchwell et al. (2006, 2007) surveys. The organizers of the MWP attribute this large increase to the 10,000-fold increase in human classifiers and the use of 24 μm data to further emphasize bubble interiors. They estimate the MWP catalog to be 95% complete, based on the falloff rate at which new objects were discovered. Still, they cautioned that this catalog is heterogenous, and probably affected by hard-to-measure selection biases.

Zoom In Zoom Out Reset image size

In terms of accuracy, humans still outperform computers in most image-based pattern recognition tasks (e.g., Zhang & Zhang 2010). Because of this, morphologically complex structures in the ISM (including supernova remnants, outflows, bubbles, H ii regions, molecular and infrared dark clouds, and planetary nebulae) are still traditionally cataloged manually. Human classification has several disadvantages, however.

First, human classification is time-consuming, and people-hours are a limited resource. Even by enlisting large communities of citizen scientists, data from next generation surveys will be too large to search exhaustively. For example, the >35,000 citizen scientists of the MWP classified roughly 45 GB of image data from Spitzer. Current and next-generation astronomical data sets are many thousands of times larger than this, suggesting tens of millions of citizen scientists would be needed for similar exhaustive searches through tera- and petabyte data sets.

Second, many scientifically important tasks are not suitable for enlisting the public. Part of the appeal of the MWP is due to the fact that the Spitzer images are beautiful, contain many bubbles, and are compelling to browse through. Searches for very rare objects, or tasks where the scientific justification is less apparent to a citizen scientist, may be less likely to entice large volunteer communities. Raddick et al. (2013) considers the motivations of citizen scientists in greater detail.

Finally, manual classification is not easily repeatable, and hard to calibrate statistically. For example, it is unknown how well the consensus opinion among citizen scientists corresponds to consensus among astronomers. The MWP catalog does not include any estimate of the probability that each object is a real bubble, as opposed to another structure in the ISM.

Automatic classifications driven by machine learning techniques nicely complement human classification. Such an approach easily scales to large data volumes and is immune to some of the factors that affect humans, like boredom and fatigue. Furthermore, because algorithmic classifications are systematic and repeatable, they are easier to interpret and statistically characterize. Despite the structural complexity of the ISM, Beaumont et al. (2011) demonstrated that automatic classification algorithms can discriminate between different ISM structures based upon morphology.

There have been a few previous attempts to detect shell-like structures in the ISM using automated techniques. One approach involves correlating observational data with templates to look for characteristic morphologies. Thilker et al. (1998), for example, used position–position–velocity templates of expanding shells to search for bubbles in H i data. This technique is also used by Mashchenko & Silich (1995) and Mashchenko & St-Louis (2002). Ehlerová & Palouš (2005, 2013) pursued a more generic approach—they looked for the cavities interior to H i shells by identifying extended regions of low emission. Daigle et al. (2003, 2007) searched for H i shells by training a neural network to classify spectra as containing outflow signatures, based on a training set of spectra from known shells. They pass spectral cubes through this classifier, and look for connected regions of outflow signatures.

It is difficult to robustly classify ISM structures using templates, due to the heterogeneity and irregularity of the ISM—simple shapes like expanding spheres are often poor approximations to the true ISM. The method pursued by Ehlerová et al. is more robust in this respect, but, as Daigle et al. (2007) noted, is only sensitive to fully enclosed bubbles. Real bubbles are often broken or arc-like. In this respect, the neural network approach that Daigle et al. pursue is an improvement, since neural networks make fewer hard-coded assumptions about the defining characteristics of a bubble, and instead use the data to discover discriminating features. There are a few drawbacks to this approach. First, Daigle et al.'s neural network classifies a single spectrum at a time, with no complementary information about the spatial structure of emission on the sky. Second, the neural networks in Daigle et al. (2007) are trained on data from 11 bubbles from the Canadian Galactic Plane Survey, and thus potentially biased toward the peculiarities of this small sample. Finally, neural networks have a large number of tunable parameters related to the topology of the network, and the optimal topology has a reputation for being difficult to optimize.

Compared to these efforts, the approach we pursue here is similar to the work of Daigle et al. (2007). However, our methodology differs in a few key respects. First, the MWP catalog provides a two orders-of-magnitude larger sample of known bubbles with which to train an automatic classifier. Second, we use the more recent random forest algorithm (Breiman 2001), which has a smaller number of tunable parameters that are easier to optimize. Finally, we focus on detecting morphological signatures of shells in Spitzer data—infrared bubbles are visually more prominent than H i shells, and thus potentially easier to classify.

There are several reasons to do build a bubble classifier based on MWP data.

• 1.  
  The automatic classifier is capable of producing quantitative reliability estimates for each bubble in the MWP catalog, potentially flagging non-bubble interlopers and leading to a cleaner catalog (Section 4).
• 2.  
  We can search for bubbles not detected by MWP citizen scientists (Section 5).
• 3.  
  We can treat this task as a case study for complex classification tasks in future data sets, where exhaustive manual classification will not be feasible (Section 6).

Brut uses the Random Forest classification algorithm (Breiman 2001) to discriminate between images of bubbles and non-bubbles. Random Forests are aggregates of a simpler learning algorithm called a decision tree. A decision tree is a data structure which classifies feature vectors by computing a series of constraints, and propagating vectors down the tree based on whether these constraints are satisfied. For example, Figure 3 shows a simple decision tree for classifying whether or not animals are cats. Here, the feature vectors are four properties of each animal observed: number of legs, tails, height, and mass. Two example feature vectors are shown, and propagate to the outlined classification nodes. In the context of Brut, each classification object is a square region in the glimpse survey. Each feature vector is a set of morphological properties extracted from one such region, and the decision tree predicts whether each region is filled by a bubble.

Decision trees are constructed using an input training set of pre-classified feature vectors. During tree construction, a quality heuristic is used to rate the tree. A few heuristics are common, which consider both the classification accuracy and the complexity of the tree itself—highly complex trees are more prone to over-fitting, and thus disfavored (Ivezić et al. 2014). We treat the choice of specific heuristic as a hyper-parameter, which we discuss below. Decision trees are constructed one node at a time, in a "greedy" fashion. That is, at each step in the learning process, a new boolean constraint is added to the tree to maximally increase the score of the quality heuristic. This process repeats until the quality heuristic reaches a local maximum.

On their own, decision trees are prone to over-fitting the training data by adding too many nodes. In a more traditional fitting context, this is analogous to adding too many terms to a polynomial fit of noisy data—in both situations, the model describes the input data very well, but does not generalize to new predictions. Random Forests were designed to overcome this problem (Breiman 2001). Random Forests are ensembles of many decision trees, built using different random subsets of the training data. The final classification of a Random Forest is simply the majority classification of the individual trees. Since over-fitting is sensitive to the exact input data used, each individual tree in a Random Forest over-fits in a different way. By averaging the classifications of each tree, the Random Forest mitigates over-fitting. Random Forests have proven effective in many machine learning contexts (Kuhn & Johnson 2013).

In addition to good performance and robustness against over-fitting, the Random Forest algorithm has several other advantageous properties: it is possible to interpret how important each element in the feature vector is to the classification (Section 2.7), the algorithm naturally ignores irrelevant elements of the feature vector, and it is conceptually and computationally simple.

It is worthwhile to compare the properties of learning algorithms like the Random Forest classifier to more traditional, "parametric" models commonly used in astronomy. Parametric models are often constructed by hand, and derived from an underlying astrophysical model—the fit to the power spectrum of the Cosmic Microwave Background is a quintessential example of such a model. The structure of a Random Forest model, on the other hand, is set by the clustering properties of the data themselves. The main advantage of the Random Forest approach is that the complexity of the model is to a large extent driven by the complexity of the training data. That is, Random Forests can model complex morphological differences between bubble and non-bubble images, which would otherwise be difficult to express in a manually constructed model of bubble structure. This is a crucial property when considering the irregularity of the ISM.

The individual "objects" that Brut classifies are square regions of the sky, and the goal of the classifier is to determine whether each region is filled by a bubble. Each field is identified by three numbers: the latitude and longitude of the center of the field, and the size of the field. We decompose the glimpse survey coverage into fields of 18 different sizes, logarithmically spaced from 002 to 095. At each size scale, we define and classify an overlapping grid of fields. Neighboring tiles in this grid are separated by one-fifth of the tile size.

As a preprocessing step, we extracted two-color postage stamps for each field (at 8 μm and 24 μm), and resampled these postage stamps to (40 × 40) pixels. Following a scheme similar to Simpson et al. (2012), these images were intensity clipped at the 1st and 97th percentiles, normalized to maximum intensity of 1, and passed through a square root transfer function. The intensity scaling tends to do a good job of emphasizing ISM structure, making bubbles more visible to the eye. Likewise, the (40 × 40) pixel resolution was chosen because it is reasonably small, yet has enough resolution that postage stamps of known bubbles are still recognizable as such by humans. Figure 4 shows four preprocessed fields toward known bubbles.⁵ The goal of preprocessing is to standardize the appearance of bubbles as much as possible, across different size scales and ambient intensities. All subsequent stages of Brut work exclusively with these images, as opposed to the unscaled pixel data.

The input to a Random Forest classifier is a numerical feature vector that describes the properties of each region. The ideal feature vector captures the differences between each class of objects, so that objects from different classes occupy different sub-regions of feature space. Designing an effective feature vector is often the most important determinant of classification performance.

The most obvious choice for a feature vector is simply the numerical value of the pixels in each preprocessed postage stamp. This turns out to be a poor choice, because extended objects like bubbles are characterized by correlations between hundreds or thousands of pixels, and any individual pixel is virtually irrelevant to the classification task. Machine learning algorithms generally perform better when individual elements of the feature vector are more important to the classification, and less dependent on the value of other feature vector elements.

Our feature vectors have been designed based on insights in the automated face detection literature (Viola & Jones 2001). The basic strategy is to encode a very large number (∼40,000) of generic statistics about image structures at various positions and scales. This strategy seems to be more effective than trying to tune a feature vector to the specific object being identified. While most elements in the feature vector will have no predictive power in the classification, the Random Forest classifier is capable of ignoring this information, and finding the elements in the feature vector with the most predictive power.

The following quantities are extracted from each postage stamp, and concatenated to form a feature vector.

• 1.  
  The wavelet coefficients from the Discrete Cosine and Daub4 wavelet transforms (Press et al. 2007). These coefficients represent the weights used in representing an image as a linear combination of basis functions. The basis functions in the Discrete Cosine transform are cosine functions, and thus encode information about power at different spatial frequencies. By contrast, the Daub4 wavelets are spatially compact, and thus better encode spatially isolated structures in an image. These wavelet-derived features resemble the features employed for face detection (Viola & Jones 2001).
• 2.  
  The image dot product⁶ of each image with 49 template images of circular rings of different size and thickness. Bubbles are morphologically similar to these templates, and tend to have larger dot products than a random intensity distribution. All ring templates are centered within the postage stamp. They span seven radius values ranging from 1–20 pixels, and seven thicknesses (2, 4, 6, 8, 10, 15, and 20 pixels).
• 3.  
  The byte size of each compressed postage stamp image—images substantially more or less smooth than a bubble compress to smaller and larger files, respectively.
• 4.  
  The DAISY features of the image (Tola et al. 2010). DAISY features are derived from measurements of local intensity gradients, and encode information about simple shapes in images; they are conceptually similar to edge detection filters.

We arrived at these specific feature vectors through experimentation, because they empirically yield effective classifications. Feature vector design is often the most subjective component of many machine learning tasks, and other features may better model the differences between bubbles and non-bubbles.

We found that classification performance improved if, in addition to extracting feature vectors at a single (ℓ, b, d) field, we also extracted features at (ℓ, b, 2 × d), (ℓ, b, d/2), and (ℓ, b + d/2, d)—here, d is the bubble diameter. This ensemble of fields gives the classifier context about the region surrounding a bubble.

Thus, in our approach, Brut extracts feature vector information from a region of the sky twice as large as the field it is classifying. This imposes a latitutde-dependent limit on the largest field Brut is able to classify. Because the glimpse coverage is limited to |b| < 1°, fields where |b| + d > 1° partially fall outside the survey area. We cannot derive proper feature vectors for these regions, and thus bubbles like this are undetectable by Brut. We return to this fact in Section 4.

Building a Random Forest requires providing a set of pre-classified feature vectors. The MWP catalog is a natural choice for building such a training set. However, as shown in Figure 5, the catalog as a whole contains many examples of objects that are not bubbles, and can compromise the training process. Thus, we manually curated a list of 468 objects in the MWP catalog which were clear examples of bubbles. We focus our attention on objects in the large bubble catalog, which contains 3744 out of 5106 total bubbles, with semimajor axes a ⩾ 135 (Simpson et al. 2012). We focus on the large bubble catalog since the web interface used to build the small catalog was different, and may have different selection biases. The largest bubble in the MWP catalog has a semi major axis a = 103. This upper limit is set by the widest-field images presented to citizen scientists—we discuss an implication of this cutoff in Section 4. The full training set of positive examples is obtained by extracting a field centered on each of these 468 clean bubbles, with a size equal to the major axis of each bubble.

Zoom In Zoom Out Reset image size

Training also requires an example set of non-bubbles. In principle, any random region of the sky that does not overlap an object in the MWP catalog is a candidate negative example—while some of these fields actually contain bubbles not in the MWP catalog, they are rare enough that they do not impact the training. However, we found that selecting negative examples at random was sub-optimal. The reason for this is that most random fields are nearly empty, and easily distinguished from fields with bubbles. Classifiers generated from such a data set incorrectly labeled nearly any field containing ISM structure as a bubble.

A better set of negative examples includes more fields containing structure (Figure 6). We built such a collection in a bootstrapped fashion. We began with a random set of negative fields, distributed uniformly in latitude and longitude, with sizes drawn from the size distribution of the positive training set. We trained a classifier with these examples and used it to scan 20,000 bubble-free regions. We then discarded half of the initial negative examples (those classified most confidently as not containing bubbles), and replaced them with a random sample of the misclassified examples. We repeated this process several times, but found that one iteration was usually sufficient to build a good set of training data.

Our final training set consists of 468 examples of bubbles identified by the Milky Way Project, and 2289 examples of non-bubble fields. We found that an approximate five-fold excess of negative examples yielded a more discriminating classifier with lower false-positive rates, compared to a training set with equal number of positive and negative examples. Likewise, using dramatically more negative examples compromised the classifier's ability to positively identify known bubbles. The number of positive examples was set by the number of clear bubbles in the MWP catalog itself (as judged by the lead author).

Instead of building a single Random Forest classifier, we trained three forests on different subsets of the sky. Each forest was trained using examples from two-thirds of the glimpse survey area, and used to classify the remaining one-third. The motivation for doing this is to minimize the chance of over-fitting, by ensuring that the regions of the sky used to train each classifier do not overlap the regions of sky used for final classification.

The training zones used to classify each forest are interleaved across all longitudes. For example, Figure 7 depicts how the 0° < ℓ < 10° region is partitioned—the shaded regions denote the portion of the sky used to train each forest, and the white region shows the zone each forest is responsible for classifying.

The Random Forest algorithm has a small number of tunable parameters.

• 1.  
  The heuristic used to measure classification quality when building a decision tree. We considered three heuristics: the Gini impurity, information gain, and gain ratio criteria. All of these provide similar measurements of how well a decision tree partitions a training data set—(Ivezić et al. 2014) discuss these criteria in greater detail.
• 2.  
  The number of individual decision trees to include in the Random Forest. We considered forests with 200, 400, 800, and 1600 trees.
• 3.  
  A stopping criterion, when building individual decision trees. If this criteria is set to one, new nodes will always be added to a tree, until the quality heuristic reaches a local maximum. If this number is set to a value c > 1, then it prevents a subtree from being added to a given node in the decision tree unless at least c training examples propagate to that node. This criteria can suppress over-fitting in individual trees, though it seems to be less important in the context of Random Forests. We considered values of 1, 2, and 4.

We explored the impact of these hyperparameter settings via cross validation. During cross validation, we split the training examples into two groups: a primary set with 154 bubble examples and 1462 non-bubble examples, and a validation set of 157 bubble and 10,000 non-bubble examples. We trained a Random Forest with a particular choice of hyper parameters, and measured the accuracy and false-positive rate of the classifier on the validation set. We found that a forest using the Information Gain heuristic, 800 trees, and c = 4 yielded the best performance. The quality heuristic and forest size are stronger determinants of classification performance than c. Classifier performance largely converged at 800 trees. While forests with 1600 trees did no worse than 800 tree forests, classification time was twice as slow. Once we converged on an optimal set of hyperparameters, we trained the final three Random Forest classifiers using the scheme discussed above.

The Random Forest algorithm measures the importance of each element in the feature vector during training. Intuitively, elements in the feature vector which appear in a large number of nodes (see Figure 3) are more important to the classification. Likewise, elements which appear closer to the root of the tree are more important, since a larger fraction of examples propagate through those nodes. More formally, feature importance is defined by the average improvement in the quality heuristic for all nodes associated with a feature, weighted by the number of training examples which propagate through those nodes. Feature importance is more meaningful in a relative sense than an absolute sense—features with higher importance scores have a greater impact on the behavior of the decision tree.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1 summarizes the importance of each category of features used in Brut—it shows the average and sum of the importances for all feature elements in each category. Both columns are normalized by the maximum importance value. The features derived from ring templates have the highest total importance—this is further evidenced by the fact that 19 out of the 20 most important single features are ring features (the 15th most important feature is a wavelet feature). While the average importance of the compression feature is slightly greater, there are far fewer compression features than ring features.

To classify a region after training, we compute the feature vector and dispatch it to one of the three Random Forests (depending on longitude, as discussed in Section 2.5). The forest produces a classification score between −1 and 1. This number is equal to the fraction of trees in the forest which predict the feature vector is a bubble, minus the fraction of trees which predict it is not. This score provides more information than a simple binary classification, as it gives some sense of the confidence of the classification. Furthermore, one can adjust the threshold that defines when an object is considered to be a bubble. When building a set of bubbles detected by Brut, increasing this threshold removes bubbles with lower classification scores. Ideally this increases the reliability of the catalog—since the remaining bubbles have higher scores and judged by the classifier to be clearer bubbles—at the cost of completeness.

One standard way to summarize the performance of a classifier is to plot the false-positive rate (fraction of negative examples incorrectly classified) versus the true positive rate (fraction of bubbles correctly classified) as the threshold varies. This is called the Receiver Operating Characteristic, and is shown for the three classifiers in Figure 8. The false-positive rate is measured by classifying ∼50,000 random negative fields, and the true positive rate is measured using the 468 known bubbles in the training set—note, however, that, given our partitioning strategy, each of these bubbles is classified by a forest that has not seen this example before. Also recall that the set of bubbles in this test set were selected because they are very well defined. Thus, Figure 8 overestimates the true positive rate when presented with more ambiguous bubbles. We compare the performance of the classifier on more representative bubbles in Section 3.

Zoom In Zoom Out Reset image size

The second forest has a slightly worse performance curve—at a true positive rate of 80%, the forest that classifies zone 2 has a false-positive rate about twice that of the other forests. In absolute terms, this false-positive rate of 0.0005 corresponds to six false positives out of the 13,000 test examples classified by this forest.

Using our automatic bubble classifier to conduct a full search for bubbles in Spitzer data is straightforward. The main task involves scanning the classifier across the glimpse survey footprint at different size scales, as described in Section 2.2. The classifier assigns a score to each field of view, as discussed in the previous section. To convert this continuous score into a binary "bubble/not-bubble" classification, we can choose a numerical threshold, and consider any score above this threshold to be a bubble detection. From Figure 8, the choice of this threshold sets the tradeoff between completeness (high true positive rate) and reliability (low false-positive rate). Because different scientific applications have different requirements on completeness and reliability, it is useful to set this threshold on a per-analysis basis. As we also discuss shortly, we can turn these confidence scores into probabilistic estimates that an expert would consider the region to be a bubble, by calibrating against expert classifications.

Note that Brut only attempts to determine the location and size of bubbles. Citizen scientists, on the other hand, were able to constrain the aspect ratio, position angle, and thickness of each object.

In practice, the classifier is insensitive to small adjustments to the size or location of the field. As a result, bubbles are usually associated with a cluster of fields that rise above the detection threshold (the white circles in Figure 9). We follow a simple procedure to merge these individual detections. Given a collection of raw detections, this procedure identifies the two nearest locations at similar scales, and discards the location assigned a lower score by the classifier. This process repeats until no two fields are close enough to discard. Two fields are close enough to merge if they are separated by less than the radius of either field, and have sizes that differ by less than 50%. The red circle in Figure 9 shows the result of merging detections.

Zoom In Zoom Out Reset image size

This practice works well for isolated bubbles, but occasionally fails for bubble complexes with multiple, overlapping bubbles. A more sophisticated clustering analysis may well perform better in this situation.

Of the 92 images in the expert survey, 45 were a random subset of objects in the MWP catalog (the remaining fields are discussed in the next section). Figure 10 shows the voting breakdown for these objects. Each row corresponds to a single image, and the red, orange, and blue bars show the fraction of experts who classified the image in each category. We call attention to several aspects of this distribution.

• 1.  
  Only 15% of these objects—7 out of 45—were unanimously classified as clear bubbles—this number increases to about 50% if all irregular and ambiguous objects are optimistically treated as bubbles.
• 2.  
  About 25% of the objects—12 out of 45—are judged by a plurality of experts to be irregular, ambiguous, or non-bubbles—that is, objects which are less likely to be the result of a young star sculpting the surrounding ISM. While some of these regions might be other interesting astrophysical objects like supernova remnants, most appear to be coincidental arc-like structure in the ISM, not convincingly sculpted by young stars.
• 3.  
  The most popular category for each object is chosen by 70% of experts on average. In other words, when using Spitzer images, determining the "true" nature of a typical object in the MWP catalog is subjective at the 30% level. Citizen scientists, by comparison, exhibit much lower levels of agreement; the average hit rate for these objects is 20%.

Is it possible to determine which of the remaining objects in the MWP catalog would likely be rejected as interlopers by experts? The hit rate of each object in the catalog might correlate with expert opinion, and be a useful filter. Likewise, the confidence score computed by Brut might identify interlopers. Figure 11 plots the citizen's hit rate and Brut's confidence score compared to the expert classification. The Y position of each point denotes the level of consensus for each object—the percentage of experts who chose the plurality category. Both metrics partially separate bubbles from the other categories. The hit rate is more effective at penalizing orange and red objects, which are confined to hit rate <0.2 in panel (a). On the other hand, the hit rate is ineffective at isolating bubbles with high expert consensus. The Brut score is more effective at identifying these high-consensus bubbles, which is apparent by the cluster of blue points in the upper right corner of panel (b). We also plot in Figure 11(c) a joint score—the sum of the normalized and mean-subtracted Brut score and hit rate. This combines the strengths of the individual metrics, and achieves the best separation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The expert reclassifications of this sample of MWP objects can be used to convert a raw score like the joint score to a calibrated probability that an expert would classify an object as a bubble, given that score. To achieve this, we perform a logistic regression against each of the three scores (hit rate, Brut score, and joint score). The logistic regression only considers whether the plurality category for each object is a bubble—the consensus information is not used. The gray curves in Figure 11 trace this curve, and show the predicted probability that an object is a bubble given the score on the X axis. We reiterate that this is not a fit through the (X, Y) cloud of points in the plot; rather, it is (informally) a fit to the fraction of circles which are blue at each X location of the plot. The title of each panel reports the quality of each fit, as defined by the log-likelihood:

$$\begin{eqnarray} \mathcal {L} &=& \sum _{\rm Bubbles} \log \left(P\left({\rm Bubble | Score}\right)\right) \nonumber\\ && + \sum _{\rm Non-Bubbles} \log \left(1 - P\left({\rm Bubble | Score}\right)\right). \end{eqnarray} \tag{ 1 }$$

Panel (d) in Figure 11 shows the ideal situation, where a scoring metric perfectly separates two classes of unanimously categorized objects. The joint score combines the information from citizen scientists and Brut, and comes closest to this ideal. While the logistic curves for the hit rate and joint score look similar, the latter curve is a better fit to the data as judged by the higher likelihood of the data under this model.

The remaining 47 regions in the expert survey were selected randomly from the full set of fields analyzed during Brut's full scan of the Spitzer images. A fully random sample from this collection is uninteresting, since the majority of fields are blank areas of the sky. Instead, these 47 images are approximately uniformly distributed in the Brut confidence score. We refer to these images as the "uniform" sample.

The expert vote distribution for these images is shown in Figure 12. There are more non-bubble fields in this sample, but the overall properties are otherwise similar. The average field is classified with a 72% consensus.

Figure 13 shows the Brut score for each field in the uniform sample, compared to the plurality category (color) and level of consensus (Y axis). The Brut score does a good job of separating the high consensus objects. Overall, however, the classes overlap more for the uniform sample than for the MWP sample. This is also expected, since all objects in the MWP sample have been previously identified by citizen scientists as potential bubble sites. The uniform sample does not benefit from this information, and is representative of the harder, "blind" classification task. As before, a logistic fit can be used to estimate, for objects identified in a blind search, the probability of being a bubble given the Brut score.

We intend to use the logistic fits in Figure 11 and 13 to convert Brut's score into a probabilistic estimate that an expert would judge an image to contain a bubble. A necessary condition for this is that the model be calibrated to the data—that P(Bubble|Score) is consistent with the empirical bubble fraction as a function of score, for the objects in the expert survey. Figure 14 confirms this to be the case—it divides the MWP and uniform samples into bins of probability (given by the logistic model), and shows the empirical bubble fraction for each bin on the Y axis. Here, bubbles are defined to be objects voted as bubbles by a plurality of experts. The 25%–75% posterior probability band is also shown, derived from the finite sample size assuming a binomial distribution (note that the empirical averages can fall outside this interval). The dotted one-to-one line is the expected performance of a calibrated model, with which these models are broadly consistent.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Bubbles are frequently studied in the context of triggered star formation. The material excavated by bubble-blowing stars might induce subsequent generations of star formation. Two main mechanisms for triggered star formation have been proposed: in the collect and collapse model (Whitworth et al. 1994; Dale et al. 2007), material gathered along bubble rims eventually passes a critical density and undergoes gravitational fragmentation and collapse. In the radiatively driven implosion model, the wind or ionization shock front collides with pre-existing but stable cores, and the resulting compression triggers collapse (Bertoldi 1989).

In any particular region, finding clear evidence for triggered star formation is difficult. The typical approach is to identify an overdensity of young stars within or on a bubble rim, and/or to look for an inverse age gradient with bubble radius (Deharveng et al. 2005; Zavagno et al. 2006; Koenig et al. 2008). Such an analysis is often confounded by small numbers of YSOs and ambiguous line-of-sight distances and stellar ages. Furthermore, it is often unclear whether spatial correlations between bubbles and young stellar objects (YSOs) imply a triggered relationship between the two, since star formation is a naturally clustered process (Lada & Lada 2003).

Such problems can partially be addressed by correlating YSOs with large bubble catalogs like the MWP catalog—this does not disambiguate correlation from causation, but it can overcome problems related to small-number statistics. Thompson et al. (2012) first applied such an analysis using the Churchwell et al. (2006) catalog, and Kendrew et al. (2012) later repeated this on the MWP catalog. This analysis computes an angular correlation function (Landy & Szalay 1993; Bradshaw et al. 2011), defined by

$$\begin{equation} w(\theta) = \frac{N_{BY} - N_{BR_Y} - N_{R_BY} + N_{R_BR_Y}}{N_{R_BR_Y}}, \end{equation} \tag{ 2 }$$

where N_αβ represents the number of pairs of objects from catalogs α and β with a separation of θ, B is a bubble catalog, Y is a YSO catalog, and R_B and R_Y are randomly distributed bubble and YSO catalogs. These random locations are chosen to preserve the approximate latitude and longitude distribution of each class of objects, but are otherwise uniformly distributed. Informally, w(θ) represents the excess likelihood of finding a YSO at a particular distance θ from a bubble, relative to what is expected from a random distribution of bubbles and YSOs. We normalize the angular offset θ by the radius of each bubble, such that w(θ) traces the excess YSOs as a function of offset in bubble radii. An ideal signature of triggered star formation, then, would be a local maximum of w(θ) at θ = 1.

We reproduce the Figure 15 of Kendrew et al. (2012) in Figure 20(a). This shows the angular correlation function between the MWP large catalog and rms catalog of YSOs and compact H ii regions (Mottram et al. 2007), as a function of normalized bubble radius. The main signal is a decaying correlation, indicative of the fact that star formation occurs in clustered environments. The prediction from triggered star formation is that there should be an additional peak at θ ∼ 1 bubble radius. Such a signal is not obvious in this figure, though Kendrew et al. (2012) report evidence for a peak for the largest bubbles in the MWP catalog. Likewise, Thompson et al. (2012) report a clearer peak when considering the expert-selected bubbles in Churchwell et al. (2006).

In Section 3, we demonstrated that roughly 30% of objects in the MWP catalog are interlopers—random ISM structures incorrectly tagged as bubbles. Furthermore, the relative false-positive rate increases toward giant H ii regions—precisely the regions where one might expect triggered star formation to be most apparent. Interlopers might significantly dilute bubble/YSO correlations in the MWP catalog. Fortunately, our bubble probability scores allow us to identify many of these interlopers, yielding a higher reliability catalog that still has three times as many bubbles as used by Thompson et al. (2012).

Figure 20(b) shows the same correlation function, for the subsamples partitioned by bubble probability. This reveals an additional excess from 0.5 < θ < 1 bubble radius, whose strength increases with bubble probability. This is similar to the trend with bubble size that Kendrew et al. (2012) reported, but the signal here is both stronger and due to a different subsample of bubbles—the size distribution of bubbles does not vary significantly between the three probability bins.

Star formation is a naturally clustered process. The curves in Figure 20 trace the excess of rms YSOs near bubbles relative to a purely random distribution of objects—they do not measure the excess density of YSOs relative to other star formation regions. The red "Control" curve in Figure 20 addresses this. This curve was obtained by repositioning each bubble to the location of a random rms YSO, and re-running the analysis. In other words, the curve shows the natural clustering of YSOs relative to each other. The left-most point of the curve is very large, which is an artifact of the fact that each bubble lines up exactly with the YSO it was repositioned to—this creates a strong overdensity at very small θ. Beyond this point, however, the curve resembles the P < 0.5 subset. The fact that the P > 0.9 bubbles fall significantly above the control curve indicates that YSOs are clustered around high-probability bubbles even more strongly than they cluster around one another on average. However, this clustering analysis cannot determine whether the overdensity around P > 0.9 bubbles is the result of bubble-triggered star formation, or rather an indication that bubbles form in particularly active star forming regions.