Assessing the Statistical Uniqueness of the Younger Dryas: A Robust Multivariate Analysis

During the last glacial period (c. 120-11 kyr BP), dramatic temperature swings, known as Dansgaard-Oeschger (DO) events, are clearly manifest in high resolution oxygen isotope records from the Greenland ice sheet. Although variability in the Atlantic Meridional Overturning Circulation (AMOC) is often invoked, a unified explanation for what caused these ‘sawtooth shaped’ climate patterns has yet to be accepted. Of particular interest is the most recent D-O shaped climate pattern that occurred from ∼14,600 to 11,500 years ago the Bølling/Allerød (BA) warm interstadial and the subsequent Younger 5 Dryas (YD) cold stadial. Unlike earlier D-O stadials, the YD is frequently considered a unique event, potentially resulting from a rerouting and/or flood of glacial meltwater into the North Atlantic or a meteorite impact. Yet, these mechanisms are less frequently considered as the cause of the earlier stadials. Using a robust multivariate outlier detection scheme a novel approach for traditional paleoclimate research we show that the pattern of climate change during the BA/YD is not statistically different from the other D-O events in the Greenland record, and that it should not necessarily be considered unique when investigating 10 the drivers of abrupt climate change. In so doing, our results thus confirm the ambiguity of the BA/YD’s trigger and present a novel statistical framework for paleoclimatic data analysis. Copyright statement. TEXT

overturning cell to undergo such large and rapid changes (Lohmann and Ditlevsen, 2019). While variations in atmospheric circulation, sea ice extent, and ice shelf formation/collapse have all been hypothesized as triggers (Li and Born, 2019), a unifying theory has yet to emerge (Lohmann and Ditlevsen, 2018). Given that D-O events provide compelling evidence that 25 the Earth's climate can rapidly switch from one state to another, it is imperative that we determine the causes of this variability if we are to accurately predict future climate.
In the original work of Dansgaard (Dansgaard, 1985), the most recent saw-tooth shaped interstadial stadial sequence of climate change since 120,000 years ago, associated with the Bølling-Allerød warming and Younger Dryas stadial (abbreviated here to BA/YD) from ∼14,600 to 11,700 years BP was labeled as D-O event 1 (Figure 1). Since then, however, a growing body 30 of geological evidence attributing the Younger Dryas cooling to a glacial outburst flood and/or a change in glacial meltwater drainage patterns to the ocean (Broecker et al., 1989;Clark et al., 2001;Keigwin et al., 2018) has sometimes led to this episode being treated as a unique event, rather than as one of the D-O stadials (Li and Born, 2019). Evidence that the YD cooling also might have coincided with a meteorite impact capable of 'blocking out' incoming solar radiation has helped bolster this notion (Firestone et al., 2007). In this paper, we use a multivariate outlier method to re-examine the extent to which the 35 BA/YD should be considered 'unique' in the context of the other D-O events. The motivation for this study derives from the remarkable similarities in shape (see further definition of shape below i.e. deviation from center over its timespan) between the BA/YD and other D-O events within the last 120,000 years, which leads us to question the uniqueness of the BA/YD in the Greenland record. Our approach is particularly novel for traditional paleoclimate research and we suggest it may be useful in future studies with similar aims argue for the increased implementation of similarly robust statistical methods in future 40 research. Indeed, the application of our statistical method to assessing the 'BA/YD uniqueness' is just one example given that it is exceptionally useful when synthesizing and compare climate records that are both uncertain and complex. It should also be noted that modeling studies suggest that forced and unforced AMOC variations have very similar signatures (Brown and Galbraith, 2016), so the outlier detection technique is not aimed to assess the qualities of D-O events as they result from specific triggers, but rather to provide a framework for situating the BA/YD within a broader context of many other D-O events, each 45 of which may (or may not) have the same underlying trigger.

Methods
To study abrupt decadal-to-multidecadal changes in climate associated with each of the Dansgaard-Oeschger events, we examined published changes in oxygen isotope ratios (δ 18 O) and methane (CH 4 ) from the NGRIP Greenland ice cores (Rasmussen et al., 2014;Baumgartner et al., 2014) and δ 18 O and carbon dioxide (CO 2 ) changes from the EDML, WAIS, Siple Dome, and 50 TALDICE ice cores recovered from Antarctica (Barbante et al., 2006;Bereiter et al., 2015) that span the last 120,000 years of Earth's climate history. Our NGRIP and EDML records both use the GICC05 age scale, and EDML was chosen because its spatial resolution and record length is comparable with the Greenland ice core records. Indeed, the snow accumulation at EDML is two to three times higher than at other deep drilling sites on the East Antarctic plateau, so higher-resolution atmosphere and climate records can be obtained for the last glacial period, making the EDML core especially suitable for studying 55 decadal-to-millennial climate variations in Antarctica. Including EDML δ 18 O allows us to observe changes in NGRIP δ 18 O as distinct in location but similar in meaning. This allows us to make conclusions about how the BA/YD may not have been a unique event in Greenland, but perhaps was so in the southern Atlantic. In general, our choice of records is based on those with the highest spatial resolution and tradition in the field of paleoclimatology of using these to study climate variability during both D-O events and the BA/YD. The high temporal resolution of the ice cores during the last glacial period makes them ideal 60 for use in our work. Both 18O records (NGRIP and EDML) provide local approximations of climate, CH 4 is a globally integrated signal indicative of hydrology in wetland regions (Brook et al., 2000) and CO 2 shows a strong correlation with Antarctic temperature on millennial scales (Bauska et al., 2021). Furthermore, we use both δ 18 O records (NGRIP and EDML) to provide local approximations of climate, whereas CH 4 and CO 2 are more indicative of global hydrology and temperature, respectively.
For the purposes of our study, the timing of each Dansgaard-Oeschger event is taken from the ages published in the IN-65 TIMATE (INTegration of Ice-core, MArine and TErrestrial records) dataset in Table 2 of Rasmussen et al. (2014). We then develop a stratigraphy that emphasizes the large-scale Dansgaard-Oeschger variability as follows: Firstly, all sub-interstadials in the INTIMATE record of Rasmussen et al. (2014) that are labeled by lowercase letters are very small, so we consider these to be part of the larger interstadial, and not unique events. For example, while Greenland Interstadial 1 (i.e. the BA interstadial) comprises of sub-events GI-1a through GI-1e, in our analysis this is simply treated as GI-1. A second set of sub-events in 70 the INTIMATE dataset are also denoted by decimals in Rasmussen et al. (2014). For example, Dansgaard-Oeschger event 2 in the INTIMATE dataset is separated into two sub-events, labeled GS 2.1/GI2.1 and GS 2.2/GI 2.2. Due to their generally high amplitude and tendency to span multi-centennial timescales, these sub-events must at least initially be considered as Dansgaard-Oeschger 'candidates', and thus require a more rigorous procedure to be dealt with. Firstly, we consider cases when two sub-events occur in succession and define a duration-based algorithm to determine whether each one should be considered 75 a separate Dansgaard-Oeschger event, both combined into one single event, or omitted from our analysis entirely. Figure 2 outlines this algorithmic consolidation process in the form of a flow chat using parameters x, y, and z.
Of the eight Dansgaard-Oeschger events in this period containing two sub-events -namely numbers 2,5,15,16,17,19,21, and 23 -our main analysis leads to the selection of stadial and interstadials found in Table 1. The selection of these events is based on using duration parameter choices: x = 300 yrs, y = 300 yrs, z = 200 yrs, which are at the shorter end of 80 what has previously been accepted as the length of a stadial or interstadial (e.g. Rasmussen et al. (2014)), but our results are not very sensitive to the chosen length. For example, columns 2-3 of Table 1 show that altering these parameters to (x, y, z) = (90, 100, 140) or (x, y, z) = (90, 100, 90) yields results that are 86-93% similar.
Taking D-O event 2 as an example, we observe that GI2.2, GS2.2, and GI2.1 span 120, 200, and 120 years respectively, and thus the algorithm in Figure 2 leads to the combination of GI2.2, GS2.2, and GI2.1 into a single interstadial, since the 85 sub-events are less than the parameter choices x = 300, y = 300, z = 200 respectively. In D-O event 5, however, GI5.2, GS5.2, and GI5.1 span 460, 1200, and 240 years, respectively, and thus under the same parameter choices, the interstadial-stadial choice algorithm in Figure 2 dictates that each sub-event should be treated as its own stadial or interstadial. Note that our final results differ minimally based on how sub-cycles are chosen.
Beyond ∼104 kyr BP, the CO 2 record contains only one data point for about every 500 years. Thus, to ensure the existence of a well-defined and complete record for our chose data, we restrict our analysis of the last glacial cycle to the period of 104-11 kyr BP containing D-O events 1-23. Of the eight D-O events containing sub-events within this period, our algorithm discards the second sub-event of four D-O events (i.e., 16.2, 17.2, 21.2, and 23.2), includes two second sub-events as distinct (i.e., 5.2 and 19.2), absorbs GI15.1 into the sub-stadials surrounding it, and absorbs GS2.2 into the sub-interstadials surrounding it (see Table 1 for the algorithm's decisions for other parameter values). This amounts to the consideration of 25 D-O events, four of 95 which are sub-events (i.e., events 5.1, 5.2, 19.1, and 19.2).
To initially examine the uniqueness of the pattern of climate change during Dansgaard-Oeschger event 1 (the BA/YD), we overlaid the NGRIP δ 18 O record of each D-O event over the BA/YD. To better visualize and compare the shape of each D-O event, we normalized the timescale that covers each D-O event and centered each record at its median. Here the term "normalizing" refers to rescaling the time axis stretching/shrinking the time over which each D-O event occurred to a consistent 100 number of years, and "centering" refers to positioning each D-O event at its median to observe the anomaly relative to the median in time space by subtracting the median value of each record during each D-O event. We then narrowed down the number of D-O events (including BA/YD) by visually selected those events that most closely resembled the pattern of NGRIP δ 18 O during the BA/YD. This process does not apply is not related to the statistical analysis that follows in which all 25 events were included, but provides preliminary evidence for the fact that the BA/YD's shape in the context of the Greenland records 105 is not unique in terms of the general shape of many D-O events.
In our second, arguably more important line of analysis, we investigate rigorously the shape (i.e., time evolving variability) of each of our chosen climate records (NGRIP δ 18 O, EDML δ 18 O, compiled Antarctic CO 2 , and NGRIP CH 4 ) during all Dansgaard-Oeschger events (25 including the BA/YD) using a robust principal component based outlier detection method entitled PCOut (Filzmoser et al., 2008). To perform this investigation, we first calculated (i) the magnitude of change from 110 interstadial to stadial (peak-to-trough analysis), (ii) the rate and direction (slope) of change of each record during each stadial, and (iii) the median value of each record during each stadial. Measurements (ii) and (iii) were considered only for stadial periods due to the following exploratory findings: When measurement (ii) was taken on the interstadial data, no significant variation from the typical sawtooth pattern was found, and thus stadials appeared to be a more ripe area of study. When measurement (iii) was taken on the interstadial data, it mirrored (iii) from the stadial periods, and was thus redundant. The peak-to-trough 115 analysis measured the amplitude of change from the interstadial to the stadial by calculating the difference between the mean of the warmest interstadial points and the mean of the coldest stadial points for each D-O event in the NGRIP δ 18 O record. To ensure that the peak interstadial warmth and maximum stadial cooling are selected, the mean values are calculated using only the upper and lower 10% of the δ 18 O values, respectively ( Figure 3). We calculate this peak-to-trough measure for the other three records by taking the difference of the mean of its values within the time window of NGRIP δ 18 O's maximum (minimum) 120 10% interstadial (stadial) values, acknowledging that some age uncertainty between the records may be present. However, the nature of these age uncertainties is not well known, so we use the aforementioned average of 10% maxima and minima as a robust protection against any age uncertainty. For some records, where no data exist in a given short interstadial (stadial) time window, we take the maximum (minimum) of a 300 year moving gaussian filter (250yr for CO 2 , in order to give higher weight to each of the sparser points in the dataset). While not ideal, this is the best approximation that our data limitations can offer. 125 We estimated the linear slope, and thus overall rate and direction of change, of each record during each of the stadials using ordinary least squares (OLS) regression. In many cases the NGRIP δ 18 O record behavior during stadial periods is generally flat, so records with highly negative peak-to-trough measurements and stadial slopes close to zero are a good indicator of Dansgaard-Oeschger event behavior. Finally, the median of each record for each stadial in our analysis was calculated. The values of each of these metrics for each record across all 25 chosen D-O events are shown in Tables 2 and 3. It should be noted 130 that all of the above stated measurements are robust to age uncertainties of at least 100 years. Thus, we can be fairly confident that age and delta age uncertainties will not wildly skew our results.
PCOut A robust principal component based outlier detection method, entitled PCOut, based on Filzmoser et al. (2008) was then applied to the results from our three metrics to test if the BA/YD is statistically different from other D-O events. This algorithm is proven to be efficient in high dimensions and especially effective in identifying location outliers, which is ideal 135 for the data considered here our data. We accept PCOut's slightly higher amount of false positives (i.e., higher size) than other algorithms on the basis that its extremely low level of false negatives (i.e., high power) is more important for this study since the areas in which the Younger Dryas is not unique is of particular interest. PCOut differs from typical principal component analysis schemes in two ways: 1) it robustly transforms the data before extracting principal components, and 2) it computes two measures of variance: one based on center location and the other based on scatter. In short, PCOut first shifts an n × p 140 data array by its variable-wise median and scales it by its variable-wise median absolute deviation (MAD), both of which are more robust (i.e., error resistant) estimators of location and scale (respectively) than sample mean and variance (Filzmoser et al., 2008). In our case, we let n = 25 correspond to the number of D-O event observations, and let p = 12 variables denote the result of obtaining the three aforementioned metrics on each of the four records (NGRIP δ 18 O, CH 4 , and δ 18 O and CO 2 from Antarctica). PCOut then performs a standard principal component analysis (PCA) procedure to the transformed data 145 that retains the first p * components contributing 99% of the data's variance, and subsequently shifts and rescales the principal components once again by their new median and MAD. For an estimate of location exceptionality, PCOut is programmed to weight each of these resulting components z * ij by the following robust measure of kurtosis, and then computes a robust Euclidian distance RD i for each of the n data points using these weights, where W = p * j=1 w j , 150 the total weight, and the z * ij are the location shifted and rescaled principal components (visualized in panel 1, Figure 4): (2) This is followed by a further transformation to acquire the final robust distances d i , where χ 2 p * ,0.5 is the 50th percentile of a chi-squared distribution with p * degrees of freedom: .
(3) 155 These d i 's represent the degree of separation each of the n data points (corresponding to D-O cycles) experiences from the center of the data, where each RD i is a robust calculation of the ith data vector's distance from its variable-wise median.
Dividing by the median of the RD i 's as in eq. 3 measures how much each RD i deviates from the median of all such distances.
To evaluate these distances as outliers or non-outliers, each data point is assigned a weight a i based on its distance d i such that higher distances receive a smaller weight so as to avoid outlier masking (visualized in panel 2, Figure 4): For this a i weight, the parameter M is the 1 3 quantile of the distances {d i }, and Finally, PCOut defines another metric b i for each data point that uses the same exact procedure as that of a i , minus the kurtosis weighting step (i.e., unweighted euclidian distance substitutes for eq. 2). Since the non-kurtosis weighted d i 165 are proven to follow χ 2 p * relatively closely, (M 2 , c 2 ) = (χ * p * ,.25 , χ * p * ,.99 ) when eq. 4 is applied to calculate b i . In the final test (visualized in panel 3, Figure 4), outliers are then defined as data points where Essentially, eq. 6 combines the result of kurtosis-weighted and non-kurtosis weighted weights a i and b i in order to provide a comprehensive measure of exceptionality. The constants 0.25 and 1.25 are added in eq. 6 to ensure that the left of side of this Further, its use of robust statistical estimators suits our study constructed dataset well in that the metrics calculated are subject to high uncertainty. For a graphical representation of PCOut's data transformations spanning eq.'s 2-6, see Figure 4.

175
PCOut is not applicable to single variable data because principal component analysis is not a valid procedure for p = 1, so outliers in this case are determined using a simpler criterion: if some point x i is less than the first quantile of the data minus MAD(x 1 , . . . , x n ) or greater than the third quantile of the data plus MAD(x 1 , . . . , x n ), it is considered an outlier.  Table 4, which indicates the subsets of records and metrics for which the BA/YD is an outlier or not -"YES"

Results
("NO") means that the BA/YD is (not) an outlier within that subset. Note we are specifically looking for subsets of records that PCOut identifies as non-outliers, and are not seeking compare the results of such subsets to one another. Thus, the temptation 200 to tally the results should be resisted -this is not a statistically sound way of determine outlier behavior on a large scale.
Beginning with single variable results, we find that all but two cells in the median column (column 4) of Table 4 exhibit outlier behavior. This result is unsurprising given that the BA/YD occurs during a period of overall warming closer to the Holocene and thus higher percentages of all chemical records compared to other D-O events, which all occurred during the coldest stretches of the past 120 kyr. Median measurements for other D-O cycles on the edges of the last glacial period also 205 harbor a proportionally higher level of median measurements due to their temporal proximity to warmer periods before and after the last glacial period. In fact, we find that the this same table of median-only PCOut results for D-O events 2, 20, and 23 harbor 47%, 60%, and 53% outliers, respectively, which indicates that we can consistently expect subsets of measurements including the median to be greater for D-O events near the beginning and end of the last glacial period. Thus, we attribute the a portion of the BA/YD's outlier behavior in variable subsets including the median to a known temperature increase during the In the single variable stadial slope column (column 3, Table 4), we find particular interest in the fact that all pairs of records including NGRIP δ 18 O (rows 6-8) do not register as outliers, while all pairs of records not including NGRIP δ 18 O (rows 9-11) do register as outliers. This exact phenomenon is also reflected in the paired peak-to-peak and stadial slope column (column 1, Our main goal is to assess the exceptionality of the BA/YD in the Greenland record. While the single record NGRIP δ 18 O section (row 12, Table 4) exhibits a mix of measured outliers and non-outliers, it must be taken into account that all variable subsets in this row that cause the BA/YD to become an outlier contain the median measurement, which, as previously stated, contributes significantly to the BA/YD's outlier behavior due to known warming leading up to the Holocene.

220
The single record NGRIP CH 4 row (row 15, namely, no outlier behavior in metrics other than the pure median exists. Since much of the ice-core based knowledge generated on the Younger Dryas relies on these two records, A lack of outlier behavior in shape is a major result, and confirms our analysis of Figure 5. Thus, the fact that the pair of behaviors of δ 18 O and CH 4 during the BA/YD is not unique compared to the pair 230 of behavior for these records in other D-O events is a stronger conclusion than if we were to restrict this analysis to only one record from Greenland. We observe this in the scatterplots of Figure 6, which plot the value of pairs of metrics for the Greenland shape measurements across all 25 D-O events, and clearly indicate that for each such pair, the BA/YD is within the natural scatter range of all other D-O events. In particular, notice that for the paired peak-to-trough scatterplots (third panel down from first column of Figure 6), the distribution of points roughly forms a ring of which BA/YD is a part, since there is  Furthermore, it should be noted that this same PCOut procedure was applied to all other D-O cycles under consideration. No clear trends were found amongst these events as a whole, but two main observations can be made. Firstly, and most importantly, it should be noted that no event exhibited non-outlying behavior isolated within the Greenland shape data as did the BA/YD.
This suggests that later research may be successful in proving that the BA/YD's data from Antarctic sources is uniquely shaped as compared to other D-O events. Secondly, we find that most data subsets for D-O events in the middle of the observed our 245 timescale (c. 49-28 ky BP) are overwhelmingly non-outlying, whereas data subsets associated with D-O events on the tails of the observed our timescale are more sporadically outlying and non-outlying. From this we might conclude that the period spanning D-O 3 to D-O 13 generally consisted of regular and "typical" D-O events, whereas D-O events not in this period either have average higher temperature (as previously discussed) or other inconsistencies. This observation does not, however, negate the conclusion that the BA/YD's Greenland data is non-outlying.

250
Note that we use the stadial/interstadial length parameters (x, y, z) = (300, 300, 200) to choose 25 D-O cycles for this section, but different parameter choices that output 28-30 D-O cycles for analysis (see Table 1) yield results that are 86-93% similar across all D-O cycles. We find this by applying PCOut to all subsets of the 28 and 30 cycle versions of our algorithm output created by implementing the parameters in columns 2 and 3 of Table 1, then comparing the results to our chosen version.

255
The aim of this study is to precisely and robustly classify the record-based qualities that would render the BA/YD a unique climate event in the context of other abrupt episodes of climate change during the last 120,000 years, known as Dansgaard-Oeschger events. If the BA/YD is to be excluded from the list of D-O events, or assigned its own particular set of triggering mechanisms there must be some statistically sound reason for doing so. Using four chemical records commonly included in assessments of general D-O behavior -δ 18 O and CH 4 from NGRIP, Greenland, δ 18 O from EDML, Antarctica, and compiled 260 CO 2 from multiple Antarctic records -we refrain from performing traditional cross correlation analysis to test for lags, and instead employ a more holistic approach that captures the shape of each D-O cycle in terms of multiple variables. Three measurements to characterize both the location (median) and shape (peak-to-trough difference, stadial slope) of each chemical record for each D-O cycle are taken, and inputted into a robust principal component analysis algorithm (PCOut) to test for outliers.

265
Our main result is as follows: the observed data for the BA/YD is not a unique compared to that of the other D-O events recorded in the Greenland ice core record, other than the fact that its median δ 18 O levels are higher due to its proximity to deglacial warming into the Holocene. The higher increase in median δ 18 O is also not unique to the BA/YD, as D-O events 2, 20, and 23 exhibit a similar phenomenon, which we attribute to their occurrence proximal to long term global climate fluctuations.
The non-uniqueness of the BA/YD's shape is clearly indicated by the statistical indistinguishability of the changes in the to explain D-O events, rather than relying on the meltwater hypothesis. Indeed, the role of meltwater forcing in triggering the YD has been questioned a number of times since it was first proposed by W. Broecker and others in 1989. For instance, the likely that freshwater forcing was less during this period (Abdul et al., 2016), making it difficult to explain how the overturning circulation remained weakened for the 1000 year duration of the YD stadial (Renssen et al., 2015). In addition, the termination of the YD, and subsequent rapid warming into the Holocene coincide with a time of increasing meltwater runoff to the North Atlantic (e.g. Fairbanks (1989)) as the Laurentide Ice sheet over North America finally collapsed.
Data availability. Metric measurements for this project can be found in Tables 2 and 3. The data for all the proxy records analyzed are 285 available in the published literature cited in the references.
Author contributions. HN Author 1 conducted all data analysis, figure creation, and manuscript drafting. AC Author 2 advised throughout the research process and edited manuscript as necessary.
Competing interests. The authors declare that they have no conflict of interest.
Vertical lines indicate the main 25 interstadial-to-stadial transitions used in this study, labeled by number from Rasmussen et al. (2014).     . PCOut's main outlier decision steps labeled by equations 2, 4, and 6, with arrows representing transitions between sequential steps of the algorithm . After completing PCA to produce centered and rescaled components z * ij , eq. 2 (visualized in the top left panel) calculates the "distance" RDi of the ith component vector from zero with sums of squares. After rescaling the RDi's to create new "distances" di, we calculate quantities ai, bi based on the function in the bottom panel (eq. 4), such that large distances di translate into smaller values of ai, bi.
Finally, the top right panel illustrates the region by which PCOut classifies ai, bi as indicative of the ith datapoint being an outlier or not (eq. 6).  . Greenland shape scatter grid. Of the 60 multivariate subsets analyzed for outliers, the above represents the essence of our results.
Using both peak-to-trough stadial slope and measurements of NGRIP δ 18 O and CH4 (left), we observe minimal outlier behavior from the BA/YD in each pair of the four variable system, which indicates that the shape of the BA/YD's Greenland records is not unique in of itself.
Histograms along the diagonal plot the corresponding single-variable distribution, where the horizontal location of the YD's measurement is in purple on a normalized scale (i.e., the height of the purple bar is 1). Table 1. Three different stadial choice situations for different duration parameter choices in Figure 2. The choice results in the first column ((x, y, z) = (300, 300, 200)) represent our preferred choices for statistical analysis. Note that basing analysis on the stratigraphic choices represented in the second or third columns yields 86-93% similarity in results.  Table 4. PCOut BA/YD Results. "YES" ("NO") indicates that the BA/YD is (not) an outlier in the subset indicated by the row/column combination in which it's located. Rows refer to the record(s) under analysis (δ 18 O= NGRIP δ 18 O, compiled Antarctic CO2 = CO2, EDML= EDML δ 18 O, and CH4 = NGRIP CH4), and columns refers to the metric(s) applied to those records (P2T=peak-to-trough, Slp= stadial slope, and Med= median). This amounts to an n × p-variate input into PCOut, where n denotes the number of records included and p denotes the metrics applied to all such records.